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Concurrent  MR-NIR  Imaging  for 
Breast  Cancer  Diagnosis 

Birsen  Yazici 


I.  INTRODUCTION 

Near  infrared  (NIR)  diffuse  optical  imaging  provides  quantitative  functional  information  from  breast 
tissue  that  can  not  be  obtained  by  conventional  radiological  methods.  NIR  techniques  can  provide  in 
vivo  measurements  of  oxygenation  and  vascularization  state,  the  uptake  and  release  of  molecular  contrast 
agents  and  chromophore  concentrations  with  high  sensitivity.  There  is  considerable  evidence  that  tumor 
growth  is  dependent  on  angiogenesis  [25]-[27],  and  that  tumor  aggressiveness  can  be  assessed  from  its 
increased  number  of  new  vessels  and  reduced  oxygenation  state  relative  to  normal  breast  tissue  and 
benign  breast  lesions  [28]-[30].  NIR  diffuse  optical  tomographic  (DOT)  methods  has  the  potential  to 
characterize  angiogenesis  related  vessel  density  as  it  measures  the  total  hemoglobin  concentration  and 
provide  the  ability  to  differentiate  between  benign  and  malignant  lesions  based  on  oxygen  saturation. 
Furthermore,  NIR  methods  are  non-ionizing,  relatively  inexpensive  and  can  be  made  portable. 

The  diagnosis  and  management  of  cancer  involves  several  stages  where  magnetic  resonance  (MR) 
plays  a  valuable  and  growing  role.  MRI  of  the  breast  is  now  a  routine  part  of  the  clinical  care  in  many 
centers  [35]-[37].  Magnetic  Resonance  imaging  (MRI)  is  indicated  in  patients  with  inconclusive  clinical 
and/or  mammographic  examinations.  Patients  that  may  benefit  include  women  with  radiographically 
dense  breasts,  and  high  risk  potential  population  [38]-[39].  MRI  possesses  less  than  10%  contrast  for 
soft  tissue  pathology  [40],  Gadolinium  (Gd)  enhanced  MRI  offers  much  better  contrast  and  is  specific 
for  tumor  vessel  imaging.  However,  the  signal  in  the  Gd-MRI  arises  from  the  larger  vessels  as  the 
contrast  agent  is  flushed  out  of  the  vascular  bed  of  the  tumor  [41].  In  comparison,  NIR  measurements 
of  absorption  have  extremely  high  contrast.  It  was  reported  that  5%  change  in  vascular  density  as 
measured  histologically  in  ductal  carcinomas  leads  to  approximately  300%  contrast  in  NIR  absorption 
coefficients  [31].  Furthermore,  there  are  studies  suggesting  that  the  kinetics  of  contrast  enhanced  optical 
spectroscopy  provides  information  about  the  cellular  spaces  [34],  On  the  other  hand,  NIR  DOT  suffers 
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from  poor  spatial  resolution  and  as  such,  it  is  unlikely  that  NIR  imaging  will  be  a  stand-alone  screening 
method  in  the  general  population.  Therefore,  we  believe  that  the  concurrent  MR  and  NIR  imaging  brings 
together  the  most  advantageous  aspects  of  the  two  imaging  modalities  (structural  and  functional).  In  the 
future,  we  envision  that  this  multimodality  imaging  approach  will  lead  to  high  resolution  hemoglobin 
tomography  and  comprehensive  quantitative  functional  tissue  characterization  to  differentiate  malignant 
and  benign  tumors. 

In  this  project,  the  clinical  studies  are  performed  using  the  novel  MR-NIR  hybrid  time-resolved 
spectroscopy  (TRS)  imager  and  fast  Indocynine  Green  (ICG)  enhanced  spectroscopic  imager  developed  by 
Dr.  Chance,  a  Co-PI  of  this  proposal,  at  the  University  of  Pennsylvania  (UPenn),  Biophysics  Department, 
Diffuse  Optical  Imaging  and  Spectroscopy  Laboratory.  The  clinical  study  has  been  approved  by  the  UPenn 
Institutional  Review  Board  under.  The  project  plans  to  leverage  the  expertise  of  Dr.  Linda  Nunes,  M.D; 
Co-PI  of  this  proposal,  at  Drexel  School  of  Medicine,  Radiologic  Sciences  Department,  in  systematic 
interpretation  of  MR  breast  architecture  and  kinetics  for  diagnosis.  The  PI  and  her  collaborators  at  Upenn 
have  been  working  on  DOT  image  reconstruction  problem  and  developed  a  number  of  approaches  for 
DOT  using  a  priori  anatomical  information  in  the  past  [42].  The  PI  also  developed  statistical  tissue 
characterization  methods  for  ultrasound  breast  cancer  diagnosis  [47]-[51].  Methods,  tools  and  results 
developed  in  these  studies  are  directly  applicable  to  the  proposed  research. 

The  central  hypothesis  of  this  project  is  that  the  concurrent  MR-NIR  diffuse  optical  tomographic 
methods  coupled  with  fast  contrast  enhanced  NIR  spectroscopic  methods  provide  fundamentally  new 
quantitative  functional  and  structural  information  for  breast  cancer  tumor  characterization  and  detection. 
This  new  information  can  be  obtained  by  novel  modeling,  analysis  and  data  fusion  methods  from  the 
tomographic,  temporal  and  cellular-based  contrast  measurements,  which  exploit  fast  imaging  techniques 
together  with  TRS  tomographic  methods.  In  this  project,  we  investigate  new  methods  for  multi-modality 
high  spatial  resolution  hemoglobin  tomography,  pharmacokinetic  modeling  of  molecular  contrast  agents 
based  on  fast  NIR  spectroscopy  and  analysis  of  structural  and  functional  information  provided  by  MR  and 
NIR  imaging  methods  for  breast  cancer  detection  based  on  receiver  operating  characteristics  methodology. 
Specific  aims  of  the  project  are  as  follows: 

•  Aim  1 :  Utilize  a  priori  anatomical  information  provided  by  MRI,  to  reconstruct  3D  high  resolution 
hemoglobin,  water  and  lipid  concentration,  and  oxygen  saturation  images  directly  from  6  wavelength 
time  resolved  optical  measurements.  Evaluate  improvements  in  image  reconstruction  between  that 
of  stand-alone  NIR  and  concurrent  MR-NIR  measurements  using  water  and  lipid  images  obtained 
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from  MRI. 

•  Aim  2:  Develop  a  compartmentalized  pharmacokinetic  modeling  of  ICG,  optical  contrast  agent,  and 
extract  quantitative  parameters  that  can  characterize  tumor  metabolism  and  angiogenesis.  Compare 
ICG  kinetics  with  the  Gadolinium,  MR  contrast  agent,  kinetics  and  biopsy  findings. 

•  Aim  3:  Evaluate  accuracy  of  breast  cancer  diagnosis  based  on  the  quantitative  functional  infor¬ 
mation  extracted  from  stand-alone  NIR  system.  This  information  includes  hemoglobin,  water  and 
lipid  concentration,  optical  scatter  power  and  oxygen  saturation  images,  and  ICG  pharmacokinetic 
parameters.  Evaluate  the  added  value  of  ICG  kinetic  parameters  in  breast  cancer  diagnosis. 

•  Aim  4:  Combine  NIR  based  breast  cancer  diagnosis  features  with  the  systematic  MR  breast  archi¬ 
tecture  and  kinetics  interpretation  model  developed  by  Dr.  Nunes,  M.D,  Co-PI  of  this  proposal,  to 
evaluate  the  sensitivity  and  specificity  of  concurrent  MR-NIR  imaging  method.  Compare  results  with 
that  of  stand-alone  MR  and  NIR  results. 

In  the  following  sections,  we  will  provide  detailed  description  of  our  current  research  in  line  with  the 
statement  of  work  (SOW)  and  the  aims  outlined  above.  For  the  period  of  June  1st,  2004  to  May  31st 
2005,  SOW  includes  only  the  first  two  aims  of  the  project. 

II.  BODY  :  Aim  1-  MR-Guided  Diffuse  Optical  Tomographic  Image  Reconstruction 

The  SOW  with  regard  to  Aim  1  includes  the  following  specific  tasks: 

•  Task  1.  Develop  a  portable  toolbox  to  process  MRI  images  to  extract  prior  information  that  will  be 
used  to  constrain  the  optical  image  reconstruction.  This  toolbox  will  perform  3D  MR  image  seg¬ 
mentation  using  standard  methods  and  extract  the  following  information:  Average  optical  properties 
of  major  breast  tissue  types,  spatially  varying  optical  grid  structure.  0-2nd  month.  Completed. 

•  Task  2.  Modify  our  existing  optical  image  reconstruction  algorithms  that  are  based  on  perturbation 
model  to  incorporate  diffusion  equation  in  an  iterative  fashion.  2-8th  month.  Completed. 

•  Task  3.  Perform  direct  reconstruction  of  hemoglobin  concentration,  oxygen  saturation,  water  and 
lipid  instead  of  optical  absorption  and  scattering  coefficients.  6- 10th  month.  Completed. 

•  Task  4.  Evaluate  improvements  in  image  reconstruction  in  terms  of  spatial  resolution  and  quantitative 
accuracy  between  the  stand-alone  NIR  and  concurrent  MR-NIR  measurements  using  the  water  and 
lipid  images  obtained  from  MR  images.  10-12th  month.  Completed. 
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A.  Aitttl  -  Taskl: 

Figure  1  shows  a  snapshot  of  the  MR  a  priori  information  processing  toolbox  that  was  developed  by 
a  team  of  two  undergraduate  students  as  a  part  of  their  senior  design  project  under  the  Pi’s  supervision. 
The  code  is  programmed  in  C.  The  program  can  load  a  stack  of  clinical  MR  images  and  allows  the 
user  to  interactively  identify  the  coregistration  fiducials  mounted  on  the  MR  compression  plates,  perform 
3D  segmentation,  interactively  identify  region  of  interest,  define  an  adaptive  mesh  structure  (2D),  and 
perform  standard  2D  DOT  image  reconstruction  (2D  linear  least-squares). 


(a)  A  snapshot  of  the  toolbox  developed  by 
the  Pi’s  undergraduate  senior  design  team 
to  perform  preprocessing  on  MR  images 
and  basic  DOT  reconstruction. 


(b)  A  slice  of  the 
segmented  MR 
breast  image. 


(c)  Adaptive  mesh 
fit  to  the  segmented 
data. 


Fig.  1.  The  anatomical  and  optical  images  are  shown  on  the  left  and  right,  respectively.  Note  the  spatial  mismatch  between 
the  two  images. 


B.  Aiml  -  Task2: 

The  outcome  of  our  research  regarding  Aim  1,  Task  2  are  reported  in  the  following  publications:  [52]- 
[55]. 

Two-level  Domain  Decomposition  Methods  (DDM)  for  Diffuse  Optical  Tomography  (DOT)  Two-level 

V 

domain  decomposition  methods  are  introduced  for  the  forward  and  the  inverse  problems  of  DOT.  The 
proposed  methods  are  TMODDM(Two-level  Multiplicative  Overlapping  Domain  Decomposition  Method) 
for  the  forward  solver  and  TMSDM(Two-level  Multiplicative  Space  Decomposition  Method).  The  con¬ 
vergence  analysis  and  the  estimation  of  computational  costs  are  given.  Numerical  implementation  shows 
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that  the  proposed  algorithm  gives  less  mean  square  error  and  more  clear  image  than  the  image  produced 
by  non-DDM.  Furthermore  DDM  approach  provides  a  framework  to  perform  region  of  interest  DOT 
image  reconstruction  based  on  the  MR  a  priori  anatomical  information. 

DOT  is  formulated  by  photon  diffusion  equation  in  the  frequency  domain.  In  the  present  approach, 
we  have  used  the  finite  element  method  using  piecewise  bilinear  element  to  compute  photon  density  in 
the  forward  problem  and  the  trust  region  method  to  compute  optical  coefficients  which  is  a  solution  of 
nonlinear  constraint  minimization  problem. 

In  DDM,  the  domain  of  interest  is  divided  by  overlapping  subdomains  as  in  Figure  2.  Usually,  the 
matrix  solver  for  N  number  of  nodes  takes  iV2  or  N 3  order  of  computations.  But  at  each  subdomain 
the  computation  becomes  ( N/d )2  or  ( N/d )3  computations  for  each  iterations.  Thus  if  we  use  d  number 
of  parallel  computers,  the  computation  is  the  number  of  iteration  times  1  /(d2)  or  l/(rl3).  Although,  we 
use  just  1  computer,  the  computation  requires  the  number  of  iteration  times  1/d  or  l/(d2). 

Two-level  method  involves  computing  at  the  fine  level  with  the  aid  of  the  coarse  level  reconstruction. 
Usually,  the  fine  level  computation  decreases  high  frequency  error,  whereas  the  coarse  level  computation 
decreases  low  frequency  components  of  the  entire  domain  error.  Using  a  two-level  method  for  the  forward 
problem,  we  can  reduce  the  number  of  iterations.  In  20  x  20  meshes,  we  get  mean  square  error  less  than 
10-5  for  just  3  iterations.  It  is  known  that  this  number  of  iteration  does  not  depend  on  the  mesh  size 
and  subdomain  size.  Using  many  subdomains,  we  can  get  a  more  speed-up  of  the  computation.  In  the 
inverse  solver,  the  convergence  is  proved  assuming  the  initial  guess  is  sufficiently  close  to  the  target 
optical  coefficients.  Since  the  coarse  level  image  approximates  the  low  frequency  of  the  target  optical 
coefficients  well,  we  have  chosen  coarse  level  correction  before  DDM  step.  Numerical  experiments  shows 
that  these  two-level  DDM  converges  more  fast  than  non-DDM  in  mean  square  error.  The  tumor  image 
was  more  clear  in  two-level  DDM  than  non-DDM. 

The  number  of  iteration  obtained  by  using  d  subdomains  is  tabulated  in  Table  1.  Here,  Nn :  The  number 
of  nodes,  Ne:  The  number  of  elements  d:  The  number  of  subdomains,  Mp\  The  number  of  iteration  for 
the  forward  solver,  p :  The  order  of  computation  with  non-DDM.  For  nonsparse  full  matrix,  p  =  3,  and 
for  usual  sparse  matrix  p  is  between  2  and  3. 

Thus  if  take  more  subdomains  and  use  more  parallel  computers,  we  get  much  computation  savings. 

In  the  proposed  method,  if  we  know  a  priori  information  about  the  approximate  location  of  tumor, 
using  a  prior  information  from  secondary  modality  such  as  MRI,  X-ray  computed  tomography  (CT)  or 
by  a  posteriori  information  from  the  coarse  level  image,  we  get  better  results  by  updating  the  coefficients 
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Method 

Computational  cost 

non-DDM 

0(Np)+0(Np) 

TMSDM  on  1  computer 

MFd}-pO{Np)  +  d}~pO{Np) 

TMSDM  on  d  computers 

MFd}~pO{Np)  +  d~vO{Nf) 

TMSDM  on  d 2  computers 

MFd'pO(Np)  +  d~pO{Np) 

TABLE  I 

THE  COMPARISON  OF  THE  COMPUTATIONAL  COSTS  FOR  THE  NON-DDM  METHOD,  TMSDM  FOR  ONE,  d ,  AND  d2 

COMPUTERS. 


(a)  (b) 

Fig.  2.  (a)  4  x  4  domain  decomposition  with  20  x  20  mesh,  (b)  3  x  1  domain  decomposition  with  24  x  12  mesh. 

the  special  region  of  interest.  Updating  optical  coefficients  only  in  fourth  subdomain  in  Figure  3(e),  we 
have  a  better  results  than  the  image  of  TMSDM  on  all  the  subdomain  in  Figure  3  (d)  or  the  image  of 
non-DDM  in  Figure  3  (f). 

1)  Numerical  results:  The  proposed  method  is  implemented  and  compared  with  non-DDM  in  the 
following  settings:  ji's  =  8cm-1,  /ja  at  background  =  0.05cm-1,  f. ia  at  tumor  =  0.2cm-1,  u  =  2n  * 
100MHz,  c  =  3  *  1010  cm/sec,  a  —  1. 

The  construction  of  diffusion  coefficient  k  when  /i„  is  fixed  is  investigated  in  [53].  And  the  recon¬ 
struction  of  na  when  /z'  is  fixed  is  implemented  in  Figures  3,  and  as  in  [54],  [55]. 

The  proposed  method  gives  better  result  than  non-DDM  as  in  figure  3(d),  as  compared  to  figure  3(f). 
And  we  also  get  the  results:  Tikhonov  parameters  with  a  =  0.01  gives  better  results  than  no  regularization 
as  in  figure  3(d)  as  compared  to  figure  3(c).  Detailed  simulation  studies  can  be  found  in  [53],  [54],  [55]. 
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(a)The  image  of  /ia  to  be  reconstructed.  White 
region  represent  background  tissue  with 
tia  =  0.05cm-1  and  black  region  represent 
anomaly  with  jia  =  0.2cm-1. 


(b)The  location  of  20  sources  and  20  detectors. 


(c)The  reconstructed  image  of  the  absorption 
coefficient  using  Algorithm  TMSDM  with  a  —  0 
resulting  in  L2  norm  error  0.0949. 


(d)The  reconstructed  image  of  the  absorption 
coefficient  using  Algorithm  TMSDM  with 
a  =  10-2  resulting  in  L2  norm  error  0.0824. 


(e)The  reconstructed  image  of  the  absorption 
coefficient  using  Algorithm  TMSDM  only  on 
fourth  subdomain  with  a  —  10-2  resulting  in  L2 


(f)The  reconstructed  image  of  the  absorption 
coefficient  not  using  Algorithm  TMSDM  with 
a  =  10-2  resulting  in  L 2  norm  error  0.0926. 


norm  error  0.0755. 

Fig.  3.  The  reconstruction  of  the  absorption  coefficient  on9  =  [0, 6]  x  [0, 6]  with  2x2  subdomains,  20  detectors,  20  sources, 
and  400  pixels. 
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C.  Aiml  -  Tasks  3  and  4: 

1)  MR-Guided  Hierarchical  Bayesian  Formulation  of  the  Inverse  DOT  Problem:  The  outcome  of  our 
research  regarding  Aim  1,  Task  3-4  are  reported  in  the  following  publications:  [2],  [4],  [42],  [43],  [44], 
[45],  [46]. 

DOT  poses  a  typical  ill-posed  inverse  problem  with  limited  number  of  measurements  and  inherently 
low  spatial  resolution.  Both  the  low  spatial  resolution  and  ill-posedness  of  the  inverse  problem  require 
the  use  of  a  priori  information  with  high  spatial  resolution.  In  this  context,  we  developed  a  hierarchical 
Bayesian  approach  to  improve  spatial  resolution  and  quantitative  accuracy  by  using  a  priori  information 
provided  by  a  secondary  high  resolution  anatomical  imaging  modality,  such  as  MR  or  X-ray.  In  such 
a  dual  imaging  approach,  while  the  correlation  between  optical  and  anatomical  images  may  be  high, 
it  is  not  perfect.  The  proposed  hierarchical  Bayesian  approach  allows  incorporation  of  partial  a  priori 
knowledge  about  the  noise  and  unknown  optical  image  models,  thereby  capturing  the  function-anatomy 
correlation  effectively.  In  other  words,  the  hierarchical  Bayesian  approach  tackles  with  the  fact  that  the 
correlation  between  the  anatomical  and  optical  images  may  be  low.  For  example,  there  may  be  regions 
in  the  optical  image  that  do  not  have  any  anatomical  counterparts  such  that  the  tumor  may  be  apparent 
in  the  optical  image,  but  may  not  have  a  corresponding  signature  in  the  anatomical  image.  Furthermore, 
average  optical  coefficients  extracted  from  anatomical  images  may  be  significantly  different  than  the  true 
optical  coefficients  of  tissue.  As  a  result,  the  assumption  of  strong  optical-anatomy  correlation  may  cause 
undesirable,  erroneous  bias  in  optical  image  reconstruction.  Therefore,  more  flexible  prior  models  are 
needed  to  properly  represent  optical-anatomy  correlation.  The  hierarchical  Bayesian  framework  affords 
such  a  flexibility  in  designing  prior  image  and  noise  models.  In  the  hierarchical  Bayesian  framework,  one 
can  formulate  the  inverse  problem  in  multiple  stages  where  each  stage  includes  information  about  the 
unknown  parameters  of  the  preceding  stage.  The  first  stage  of  the  hierarchy  includes  the  data  likelihood 
and  the  first  stage  of  the  image  prior,  which  are  comprised  of  statistical  models  for  the  noise  and  optical 
image,  respectively.  These  models  include  parameters  associated  with  noise  and  image  statistics,  which 
are  not  known  precisely  in  practice.  These  unknown  parameters  are  referred  to  as  hyperparameters,  which 
can  be  regarded  as  random  variables.  The  succeeding  stage  of  the  hierarchical  formulation  incorporates 
a  priori  information  about  the  hyperparameters  in  the  form  of  prior  distributions  -  so  called  hyperpriors 
-  defined  on  the  hyperparameters.  The  incorporation  of  the  second  stage  concludes  the  design  of  the 
two-level  hierarchical  noise  and  image  models. 

In  order  to  estimate  the  hyperparameters,  we  adapt  the  linear  conjugate  gradient  (CG)  algorithm  to 


12 


include  a  hyperparameter  estimation  step  followed  by  an  image  update.  In  this  context,  we  apply  an 
iterative  empirical  Bayesian  approach  to  estimate  the  hyperparameters,  which  in  turn  gives  the  maximum 
a  posteriori  (MAP)  estimates  of  the  hyperparameters  at  each  CG  iteration  prior  to  the  image  update. 
Hence,  the  noise  and  image  models  are  accommodated  at  each  update  of  the  hyperparameters  along  with 
the  solution  process.  We  refer  to  the  following  paper  for  further  details  [4]  on  hierarchical  Bayesian 
formulation  and  estimation  of  hyperparameters. 

In  the  following  section,  we  present  the  optical  image  reconstruction  results  based  on  the  developed 
hierarchical  Bayesian  approach  in  comparison  to  the  maximum  likelihood  estimate  of  the  same  image. 
We  consider  two  cases  to  show  the  effectiveness  of  the  hierarchical  Bayesian  approach: 

1)  The  tumor  is  present  both  anatomically  and  optically;  thus  perfect  correlation  exists  between  the 
anatomical  and  optical  images. 

2)  The  tumor  is  present  optically,  however  it  has  no  corresponding  signature  in  the  anatomical  image. 
We  show  that  the  hierarchical  Bayesian  approach  and  the  associated  hyperparameter  estimation 
effectively  captures  the  inconsistency  between  the  anatomical  and  optical  images,  thereby  providing 
reliable  optical  image  estimates. 

2D  Experiment  with  MR-simulated  data: 

We  used  the  Tl-weighted  MR  breast  image  from  [1]  to  design  a  realistic  optical  breast  model  (fig¬ 
ure  4).  The  MR  breast  image  was  segmented  into  parenchyma  and  adipose  layers  by  applying  a  simple 


Fig.  4.  The  original  MR  breast  image  with  an  artificial  tumor  inserted. 

thresholding  algorithm  with  respect  to  the  MR  image  intensity  values.  Next,  a  tumor  corresponding  to 
an  infiltrating  ductal  carcinoma  revealed  by  Gd-DTPA  (Gadolinium-diethylenetriamine  pentaacetic  acid) 
enhancement  was  inserted  (shown  in  figure  4  as  well).  Each  sub-region  was  assigned  an  absorption  value 
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as  indicated  in  [1]  (pidipose  =  0.03cm  1 ,  ^renchvma  =  0.06  cm  1  ,  ^“mor  =  0.09  cm  J)  to  obtain  an 
initial  template  (figure  5(b)).  To  simulate  a  corresponding  optical  image,  zero  mean  Gaussian  noise  was 


(a)  (b) 

Fig.  5.  The  anatomical  and  optical  images  are  shown  on  the  left  and  right,  respectively.  Note  the  spatial  mismatch  between 
the  two  images. 

added  prior  to  filtering  the  image  by  a  low-pass  filter.  The  resulting  optical  image  is  shown  in  figure  5(b). 
Note  the  quantitative  and  spatial  mismatch  along  the  boundaries  and  especially  within  the  tumor.  The 
homogeneous  diffusion  coefficient  of  the  medium  was  set  to  0.042  cm.  9  frequencies  ranging  from  0 
to  244  MHz  were  used  to  obtain  729  measurements  with  9  sources  and  9  detectors  positioned  along 
x-axis  on  opposite  sides.  The  optical  medium  was  uniformly  discretized  into  90  pixels  along  x-axis  and 
60  pixels  along  y-axis  leading  to  a  total  of  5400  lxl  cm2  pixels. 

We  performed  two  types  of  experiments  to  test  the  performance  of  the  proposed  hierarchical  Bayesian 
approach  for  this  problem: 

(i)  Tumor  present  both  anatomically  and  optically:  In  this  experiment,  the  template  extracted  from  the 
anatomical  image  shown  in  figure  5(a)  was  used  to  design  the  hierarchical  image  prior.  As  a  result,  the 
optical  image  was  segmented  into  three  sub-images  each  of  which  corresponded  to  the  labeled  images 
in  the  anatomical  image  as  shown  in  figure  5(a).  In  the  design  of  the  hyperprior  defined  on  the  mean 
(i.e.  p(/^|C)),  values  that  significantly  different  than  the  actual  mean  of  the  sub-images  were  used.  Thus, 
this  experiment  evaluates  the  robustness  of  the  proposed  method  when  the  true  statistics  of  the  optical 
image  are  significantly  different  than  the  statistics  extracted  from  the  prior  anatomical  image. 

The  reconstructed  image  and  the  sub-image  zoomed  into  the  tumor  region  are  shown  in  figures  6(a) 
and  6(c),  respectively.  For  comparison,  the  maximum  likelihood  (ML)  estimate  of  the  image  are  shown 
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in  figures  6(b)  and  6(d).  The  simulation  results  show  that  hierarchical  Bayesian  approach  leads  to 
qualitatively  better  results  and  resolves  the  tumor  more  accurately. 


(c)  (d) 

Fig.  6.  The  hierarchical  Bayesian  reconstruction  of  the  optical  image  (a)  using  the  anatomical  template  shown  in  figure  5(b) 
for  the  design  of  the  hierarchical  image  model.  Figure  (c)  shows  the  image  that  zooms  into  the  tumor  region  in  the  optical 
image  shown  in  (a).  The  maximum  likelihood  (ML)  estimate  of  the  entire  image  and  the  sub-image  focusing  the  tumor  region 
are  shown  in  (b)  and  (d),  respectively.  The  rectangular  box  in  the  figures  shows  the  actual  location  of  the  tumor. 

In  figure  7,  the  estimates  of  the  hyperparameters  associated  with  the  noise  and  image  models  are  given 
as  a  function  of  the  iteration  number.  Note  that  the  mean  value  estimates  for  each  sub-image  converge  to 
actual  values,  even  though  the  corresponding  assigned  hyperparameters  regarding  the  mean  value  deviate 
from  the  true  average  optical  values  by  at  least  15%  (see  table  II).  The  experiment  also  demonstrates  that 
the  initialization  of  the  hyperparameters  does  not  have  any  effect  on  the  performance  of  the  estimation 
(figure  7). 

(ii)  Tumor  present  optically  but  not  anatomically:  In  this  experiment,  we  removed  the  tumor  region 
from  the  template  extracted  from  the  prior  anatomical  image,  but  kept  it  in  the  optical  image  as  shown  in 


15 


0.0151 

0.015 

0.015 

I 

0.015 

II 

0.0196 

0.0196 

0.0196 

1i 

0.0345 

0.034 

0.0335 


(b)  The  estimated  standard  deviation  cn 
for  each  sub-image  vs  iteration  number  are 
shown  on  the  left.  The  estimates  for  the 
parenchyma,  adipose  and  tumor  sub-images 
are  given  at  the  top,  middle  and  bottom, 
respectively.  The  noise  scale  estimate  A  vs 
iteration  number  is  shown  on  the  right. 

Fig.  7.  The  mean  value  and  the  standard  deviation  estimates  for  each  sub-image  (sub-images  are  determined  by  the  anatomical 
image)  and  noise  scale  estimate  versus  iteration  number  are  shown.  The  thick  solid  line  shows  the  estimated  values.  The  constant 
solid  line  in  (a)  shows  the  actual  mean  value  and  the  constant  dashed  line  in  (a)  shows  the  assigned  mean  value  (fa  +  pa o) 
used  in  the  design  of  the  hyperprior  defined  on  the  sub-image  means.  The  actual  mean  values  in  these  sub-regions  are  0.032, 
0.058  and  0.076,  respectively. 

figures  8(a)  and  8(b).  As  a  result,  the  optical  image  was  segmented  into  two-sub  images.  The  objective 
of  this  experiment  is  to  evaluate  how  well  the  proposed  method  reconstructs  optical  tumors  when  they 
are  not  anatomically  present. 

The  reconstructed  images  for  this  experiment  are  given  in  figures  9(a)  and  9(c),  respectively.  The  ML 
estimate  of  the  image  is  given  in  figures  9(b)  and  9(d).  We  observe  that,  even  though  there  is  a  significant 
mismatch  between  the  optical  image  and  the  anatomical  counterpart  in  the  tumor  region,  the  hierarchical 
Bayesian  formulation  leads  to  a  qualitatively  better  reconstruction  than  the  ML  approach,  even  around 
the  tumor.  Furthermore,  the  tumor  is  better  localized  as  compared  to  the  ML  solution  and  is  not  biased 
towards  the  a  priori  anatomical  image.  The  error  in  the  localization  of  the  tumor  can  be  attributed  to  the 
source-detector  geometry.  The  propagation  of  light  along  y-direction  results  in  a  smoothing  effect  on  the 
optical  image  along  y-direction.  This  effect  is  enhanced  near  source  and  detectors  due  to  the  behavior 
of  the  solution  of  the  diffusion  equation.  The  vertically  smoothing  effect  can  be  observed  in  the  ML 
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(a)  The  estimated  mean  values  fii  for  each 
sub-image  vs  iteration  number  are  shown. 
The  estimates  for  the  parenchyma,  adipose 
and  tumor  sub-images  are  given  at  the  top, 
middle  and  bottom,  respectively. 
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(a)  (b) 

Fig.  8.  The  anatomical  template  (a)  and  the  original  optical  image  (b).  Note  that  the  tumor  is  not  anatomically  present. 

estimate  of  the  image  more  apparently  (figures  9(b)  and  9(d)).  The  smoothing  effect  can  be  suppressed 
by  incorporation  of  a  priori  information  for  the  tumor  as  in  case  (i),  where  the  tumor  is  better  resolved 
(figure  6(c)).  Further  improvement  can  be  achieved  by  employing  sources  and  detectors  positioned  along 
y-axis  as  well  as  along  x-axis. 

In  figure  10,  the  average  absorption  value  of  each  reconstructed  sub- image  vs  iteration  number  is  given. 
The  sub-images  correspond  to  those  as  indicated  by  the  actual  optical  image  shown  in  figure  5.  Note 
that  the  mean  value  of  the  reconstructed  image  in  the  tumor  region  converges  to  the  actual  value  even 
though  the  anatomical  image  asserts  that  no  tumor  exists. 

The  set  of  parameters  used  in  the  design  of  hyperpriors  for  these  experiments  and  the  actual  mean  of 
absorption  values  for  each  sub-image  are  shown  in  table  II. 

TABLE  II 

The  actual  mean  of  the  absorption  values  in  each  sub-image  and  the  parameter  set  used  in  the  inverse 
PROBLEM  FORMULATIONS  FOR  THE  MR-SIMULATED  EXPERIMENTS  I  AND  II.  /la 0  =  0.0439  FOR  THIS  EXPERIMENT.  N/A 

STANDS  FOR  “NOT  APPLICABLE”. 


The  first 

experiment 

(pi  4~  paOj$i)« 

Parenchyma 

(0.038,  0.23) 

(0.015,  0.228) 

Sub-images 

Adipose 

(0.05,  0.3) 

(0.02,  0.3) 

Tumor 

(0.09,  0.54) 

(0.036,  0.54) 

The  second 

(£i  +  HaO  ,$i): 

(0.03,  0.18) 

(0.06,  0.36) 

N/A 

experiment 

(/Lri.Ti): 

(0.012,  0.18) 

(0.024,  0.36) 

N/A 

(Mr‘“aI+Mao): 

0.032 

0.058 

0.076 
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(C)  (d) 

Fig.  9.  The  hierarchical  Bayesian  reconstruction  of  the  optical  image  (a)  using  the  anatomical  template  shown  in  figure  8  for 
the  design  of  the  hierarchical  image  model.  Figure  (c)  shows  the  image  that  zooms  into  the  tumor  region  in  the  optical  image 
shown  in  (a).  The  ML  estimate  of  the  entire  image  and  the  sub-image  focusing  the  tumor  region  are  shown  in  (b)  and  (d), 
respectively.  The  rectangular  box  in  the  figures  shows  the  actual  location  of  the  tumor. 

2)  MR-Guided  Direct  Recovery  of  Blood  Volume  and  Saturation:  A  detailed  discussion  of  the  research 
described  in  this  section  can  be  found  in  [2], 

The  estimation  of  the  absorption  coefficient  at  several  wavelengths  enables  the  provision  of  spatial  maps 
of  the  targeted  chromophores.  In  the  case  of  the  breast,  the  four  chromophores  of  potential  diagnostic 
interest  are  the  oxy-  and  deoxy-haemoglobin,  the  water  and  the  lipid.  The  concentrations  of  these  breast 
constituents  are  linearly  related  to  the  absorption  values. 

Solving  the  linear  system  on  a  voxel  basis  for  a  spectral  set  of  absorption  distributions  gives  the 
functional  maps  required  for  diagnostic  purposes,  which  can  be  termed  as  “indirect  imaging”.  However, 
indirect  imaging  requires  consecutive  solutions  of  two  linear  systems,  both  of  which  are  ill-posed.  In  order 
to  overcome  the  burden  and  minimize  the  systematic  errors  due  to  the  ill-conditioning  of  both  inverse 
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(a)  The  average  value  ( jia  in  each  recon¬ 
structed  sub-image  vs  iteration  number  are 
shown.  The  estimates  for  the  parenchyma, 
adipose  and  tumor  sub-images  are  given  at 
the  top,  middle  and  bottom,  respectively. 
The  domains  of  the  sub-images  correspond 
to  the  domains  in  the  optical  image  shown 
in  figure  8(a). 


(b)  The  estimated  standard  deviation  bi 
for  each  sub-image  vs  iteration  number  are 
shown  on  the  left.  The  estimates  for  the 
parenchyma,  adipose  sub-images  are  given 
at  the  top  and  bottom,  respectively.  The 
sub-images  correspond  to  the  sub-images 
defined  by  the  anatomical  template  shown 
in  figure  8(b).  The  noise  scale  estimate  A 
vs  iteration  number  is  shown  on  the  right. 


Fig.  10.  The  mean  values  of  the  reconstructed  sub-images  vs  iteration  number  (a).  The  estimated  values  of  the  standard 
deviation  of  the  two  sub-images  and  the  noise  scale  A  are  shown  in  (b). 


problems,  a  method  aiming  at  directly  imaging  the  functional  parameters  has  been  proposed  [32].  This 
method  takes  advantage  of  the  linear  relationship  and  formulates  the  inverse  problem  in  terms  unknown 
chromophore  concentrations  [2].  This  new  linear  system  is  also  poorly  conditioned  and  great  care  should 
be  taken  during  the  pre-conditioning  of  this  sensitivity  matrix.  In  our  case,  we  use  an  average  column 
scheme  [33].  This  specific  preconditioning  scheme  led  to  the  most  accurate  reconstructions. 

The  functional  parameters  that  we  are  interested  in  are  the  blood  volume: 

[BV\  =  [Hb]  +  [Hb02]  (1) 


and 


[Sa02]  = 


[Hb02] 


[Hb02] 


(2) 


[Hb]  +  [Hb02]  [. BV }  ‘ 

Simulation  Experiment  Measurements  were  obtained  by  solving  the  frequency-domain  diffusion  equation 
with  a  finite  difference  approach  (FDM).  We  restricted  our  simulations  to  a  two-dimensional  (2D) 
geometry  for  computational  efficiency  and  to  a  continuous  wave  data  set  type.  The  slab  thickness  was 
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6  cm  simulating  a  soft  compressed  breast.  We  placed  nine  sources  on  one  side  of  the  slab  and  nine 
detectors  on  the  other  side,  both  evenly  stretched  along  8  cm  (see  figure  11).  Two  square  inclusions  of 
1  cm2  were  simulated  in  this  model. 
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Fig.  11.  Optical  model  used  for  the  simulations  (left).The  optical  properties  displayed  here  correspond  to  A  =  750 nm. 
Functional  parameters  used  to  define  the  optical  properties  to  generate  the  synthetic  measurements  and  optical  parameters  for 
the  spectral  set  investigated.  No  [H^O]  and  [Li]  constituents  were  used  in  these  simulations  for  simplicity. 


Fig.  12.  (a)  Blood  volume  (mM)  reconstructions  and  (b)  saturation  estimates  (%)  in  the  case  of  classical  indirect  imaging. 

First,  the  estimation  of  the  functional  parameters  simulated  is  accomplished  with  a  conjugate  gradient 
descent  (CGD)  algorithm  in  the  case  of  indirect  imaging  without  using  any  kind  of  prior.  The  results  are 
displayed  in  figure  12.  Then  reconstructions  using  the  indirect  imaging  approach  within  the  hierarchical 
Bayesian  framework  [2],  [4]  are  provided  in  figure  13.  First,  an  estimation  of  the  absorption  coefficients  at 
the  six  wavelengths  was  performed  using  spatial  a  priori  information  and  a  conjugate  gradient  algorithm 
with  the  PolakRibiere  method  [4].  The  assigned  mean  values  for  the  absorption  coefficients  were  the 


Fig.  14.  (a)  Blood  volume  (mM)  reconstructions  and  (b)  saturation  estimates  (Bayesian  direct  imaging.  The  assigned  mean 

values  correspond  to  the  true  concentration  values  with  a  30%  level  of  confidence. 

simulated  ones  and  30%  of  the  assigned  mean  values  were  used  as  the  standard  deviation.  Then  classical 
spectroscopy  was  performed  on  the  resulting  optical  maps  to  solve  for  [Hb02]  and  [Hb]  concentrations. 
In  the  second  case,  the  results  using  physiological  and  spatial  a  priori  information  are  depicted  in 
figure  14.  The  assigned  mean  values  for  the  chromophore  concentrations  were  the  exact  ones  with  standard 
deviations  30%  of  the  mean  values.  The  functional  quantitative  values  retrieved  from  each  case  and  for 
the  relevant  region  of  interest  are  summarized  in  table  15. 

In  this  work  [2],  we  reported  our  first  step  towards  incorporating  physiological  and  spatial  a  priori 
information  derived  from  MRI  to  assist  DOT.  Better  estimates  of  the  main  functional  parameters  that  are 
[BV]  and  [Sa02]  were  achieved. 
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Fig.  15.  Functional  parameters  recovered  with  the  three  approaches  and  for  the  three  different  functional  areas.  The  values 
proposed  here  correspond  to  the  mean  value  of  the  entire  ROI  defined  as  the  a  priori  spatial  masks.  [BV]  in  mM  and  [Sa02] 
in  %. 


III.  BODY  :  Aim  2  -  PHARMACOKINETIC  MODELING  OF  ICG 
Aim  2  of  the  SOW  involves  the  following  specific  tasks: 

•  Task  1.  Investigate  compartmentalized  pharmacokinetic  models  for  Gd  and  ICG  to  model  flow 
dynamics  of  ICG  in  extravascular  extracellular  space  and  plasma.  0-6th  months.  Completed. 

•  Task  2.  Develop  mathematical  methods  to  extract  parameters  from  ICG  pharmacokinetic  model  to 
capture  contrast  agent  flow  kinetics.  6-9th  months.  Completed. 

•  Task  3.  Perform  simulation  and  in  vivo  clinical  experiments  to  validate  pharmacokinetic  model. 
9- 15th  months.  Ongoing. 

•  Task  4.  Compare  ICG  flow  dynamics  with  the  Gd  flow  dynamics  of  the  breast  MRI  data.  15-1 8th 
months.  Ongoing. 

A.  Aim  2  -  Task  1  and  2: 

The  outcome  of  our  research  regarding  Aim  2,  Task  1  and  2  are  reported  in  the  following  publications 
and  presentations:  [5],  [6],  [7],  [8],  [9]. 

In  order  to  utilize  a  compartmentalized  model  for  the  ICG  optical  contrast  agent  dynamics,  to  extract 
parameters  that  can  characterize  tumor  physiology,  metabolism,  and  angiogenesis  using  fast  temporal 
NIR  spectroscopy  we  developed  three  different  compartmental  models  namely,  the  four,  three  and  two- 
compartment  models,  to  model  the  pharmacokinetics  of  ICG  in  cancerous  tumors. 

The  four-compartment  model  includes  capillary  region,  interstitial  fluid  region  (ISF),  parenchymal  cell 
region  and  intracellular  binding  site  as  compartments.  Figure  16  illustrates  the  capillary  and  extracapillary 
space  relevant  to  the  four  compartment  model.The  ICG,  injected  intravenously  into  the  subject,  can  pass 
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through  from  capillary  into  reversible  binding  site  inside  the  cell  through  the  interstitial  fluid  region  and 
the  parenchymal  cell  region  [11],  [12],  [13].  Moreover,  in  advanced  tumor  stages,  the  leakiness  around 
the  tumor  vessels  is  expected  to  increase,  resulting  in  higher  permeability  rates  during  the  transportation 
of  ICG  into  the  compartments. 


Fig.  16.  A  simple  illustration  of  the  capillary  extracapillary  structure. 

In  the  three-compartment  model,  the  parenchymal  cell  and  intracellular  binding  site  compartments  are 
combined  to  form  a  single  compartment  called  parenchymal  cell.  This  amounts  to  the  assumptions  that  the 
transition  of  ICG  into  the  intracellular  binding  site  is  negligible  as  compared  to  the  other  compartments 
and  therefore  omitted  from  the  model. 

In  the  two-compartment  model,  the  tumor  region  is  assumed  to  be  composed  of  two  compartments; 
namely,  the  capillary  region  and  ISF  [10],  [14],  [15].  The  transition  of  the  ICG  to  the  third  and  fourth 
compartments  are  assumed  to  be  negligible.  Therefore  the  last  two  compartments  in  the  four  compartment 
model  is  omitted.  We  consider  transcapillary  leakage  to  occur  only  at  the  tumor  site.  We  also  assume 
that  a  small  perturbation  of  the  global  plasma  concentration  does  not  affect  the  bulk  removal.  The  block 
diagrams  of  the  compartmental  transport  and  chemical  model  of  ICG  delivery  for  four,  three,  and  two- 
compartment  models  are  shown  in  figure  17  (a),  (b),  and  (c),  respectively. 

The  set  of  differential  equations  representing  the  ICG  transition  between  the  four  compartments  are 
given  as  follows: 

The  leakage  into  and  the  drainage  out  of  plasma: 

=  kbCi{t)  -  kaCp{t)  -  koutCp(t).  (3) 

The  leakage  into  and  the  drainage  out  of  ISF: 

=  kaCp(t)  -  kbCi{t)  -  kcCi{t)  +  kdCpc{t). 


(4) 
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(a) 

Fig.  17.  Block  diagrams  of  the  (a)  four-compartment,  (b)  three-compartment,  (c)  two-compartment  models  for  the  ICG 
pharmacokinetics. 

The  leakage  into  and  the  drainage  out  of  parenchymal  cell: 

=  kcCi(t)  -  kdCpc(t)  -  keCpc(t)  +  kfCb(t).  (5) 

The  leakage  into  and  the  drainage  out  of  intracellular  binding  site: 

^ l  =  keCpc(t)-kfCb(t ).  (6) 

Physiologically,  the  equilibrium  constants  are  defined  by  the  permeability  surface  area  products  given 
as  PSp,  where  P  is  the  capillary  permeability  constant,  S  is  the  capillary  surface  area,  and  p  is  the 
tissue  density.  kout  is  proportional  to  the  flow  rate  into  and  out  of  the  capillary  and  ka  kb,  kc  ,  kd,  ke, 
and  kf  represent  intra-tissue  physiologic  effects  during  ICG  delivery  from  capillary  to  binding  site. 

The  actual  bulk  ICG  concentration  in  the  tissue  measured  by  NIR  spectroscopy,  m(t),  is  a  linear 
combination  of  the  ICG  concentrations  in  four  different  compartments. 

m{t)  =  vpCp(t )  +  ViCi{t)  +  vpcCpc(t)  +  vbCb(t),  (7) 

where  vp,  Vi,  vpc,  vb,  are  volume  fractions  of  plasma,  ISF,  parenchymal  cell  region  and  intracellular 
binding  site,  respectively. 

The  differential  equations  describing  the  transitions  between  compartments  for  three  and  two-compartments 
models  are  similar  to  the  four-compartment  model  and  described  in  detail  in  [5]. 
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After  compartmentalizing  the  region  around  the  tumor  we  developed  a  mathematical  estimation  method 
to  extract  ICG  kinetic  rates.  Using  these  rates  one  can  differentiate  between  malignant  and  benign  tumors 
and/or  different  tumor  stages. 

1)  Kalman  Filtring  for  the  Estimation  of  Contrast  Agent  flow  Kinetics:  We  introduce  a  systematic  and 
robust  approach  to  model  and  analyze  ICG  pharmacokinetics  based  on  the  extended  Kalman  filtering 
(EKF)  framework.  Kalman  filter  (KF)  is  an  optimal  recursive  modeling  and  estimation  method  with 
numerous  advantages  in  ICG  pharmacokinetic  modeling.  These  include: 

•  Effective  modeling  of  multiple  compartments,  and  multiple  measurement  systems  governed  by 
coupled  ordinary  differential  equations,  in  the  presence  of  measurement  noise  and  uncertainties 
in  the  compartmental  model  dynamics. 

•  Simultaneous  estimation  of  pharmacokinetic  model  parameters  and  ICG  concentrations  in  each 
compartment,  which  is  not  accessible  in  vivo  by  means  of  NIR  techniques. 

•  Recursive  estimation  of  time-varying  pharmacokinetic  model  parameters. 

•  Statistical  validation  of  estimated  concentrations  and  error  bounds  on  the  pharmacokinetic  parameter 
estimates. 

•  Incorporation  of  available  a  priori  information  about  the  initial  conditions  of  the  permeability  rates 
into  the  estimation  procedure. 

•  Potential  real-time  monitoring  of  ICG  pharmacokinetic  parameters  and  ICG  concentrations  in  dif¬ 
ferent  compartments  due  to  the  recursive  nature  of  EKF  estimation  method. 

The  details  of  the  extended  Kalman  filter  for  simultaneous  estimation  of  pharmacokinetic  model 
parameters  and  ICG  concentrations  in  each  compartment  is  explained  in  [5],  [6]. 

For  the  optimal  compartmental  model  order  selection  we  adopted  Bayesian  information  criteria  (BIC). 
BIC  is  a  well  known  information  theoretic  criteria,  in  which  the  optimal  model  order  is  selected  by 
minimizing  a  cost  function  to  avoid  overfitting.  The  cost  function  depends  on  the  number  of  observations, 
number  of  unknown  parameters  to  be  estimated  and  the  maximum  likelihood  function.  A  detailed  discus¬ 
sion  on  BIC  can  be  found  in  [16],  [17],  [18].  In  order  to  calculate  the  BIC  for  different  compartmental 
models,  we  first  derived  a  likelihood  function  for  the  extended  Kalman  filter.  The  derivation  is  based 
on  the  maximum  likelihood  estimation  of  the  parameters  in  the  Kalman  filtering  framework  given  as  in 
[19],  [20].  We  then  modified  this  likelihood  function  for  the  extended  Kalman  filter  estimator  for  the 
joint  estimation  of  compartmental  model  parameters  and  concentrations.  The  details  of  the  derivation  is 
provided  in  [5], 
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Fig.  18.  ICG  concentrations  measured  in  tissue  for  four  different  rats. 

B.  Aim  2  -  Task  3: 

The  outcome  of  our  research  regarding  Aim  2,  Task  3  are  reported  in  the  following  publications  and 
presentation:  [5],  [6],  [7],  [8],  [9]. 

1)  Fischer  Rat  Data  Analysis:  In  order  to  validate  the  pharmacokinetic  modeling  of  ICG  kinetics 
we  first  used  the  ICG  concentration  data  obtained  from  four  Fischer  rats  with  adenocarcinoma  tumors. 
R3230ac  adenocarcinoma  cells  were  injected  below  the  skin  into  four  Fischer  rats  3  weeks  prior  to 
measurements.  The  tumor  size  for  the  rats  ranges  in  diameter  from  5  to  30  mm.  The  ICG  concentration 
data  was  collected  with  an  MR-NIR  imager.  The  configuration  of  the  apparatus,  the  data  collection 
procedure,  and  the  details  of  the  experimental  approach  have  been  reported  in  [10],  [21],  [22].  . 

Figure  18  presents  the  ICG  concentrations  (/ xM )  from  four  different  rats.  Tumors  in  Rat  1  and  2 
are  classified  as  necrotic  because  of  their  low  tissue  oxy-hemoglobin,  low  total  hemoglobin,  and  low 
gadolinium-diethylene-triamine  penta-acetic  acid  (Gd-DTPA)  enhancement  levels.  Tumors  in  Rat  3  and 
4  are  classified  as  edematous  due  to  their  high  water  content.  It  can  be  observed  from  figure  18  that  the 
necrotic  cases  display  low  peak  ICG  concentration  values  and  slowly  rising  slopes  unlike  the  edematous 
cases  with  high  peak  values  and  sharp  rising  slopes. 

We  estimated  the  pharmacokinetic  rates  for  the  four-,  three-  and  two-compartment  models.  The  results 
are  given  in  tables  III,  IV,  and  V,  respectively.  The  estimated  pharmacokinetic  rates  for  all  compartmental 


26 


TABLE  III 


Four-Compartment  Model:Estimated  pharmacokinetic  parameters  using  EKF  algorithm 


t(4) 

roa 

^4) 

1.(4)  1.(4) 

roc  ti,d 

*.(4)  1.(4)  k(4) 

fte  Kf  * 'out 

(sec-1 10-2) 

(sec-110-2) 

(sec-1 10-2)  (sec-1 10-2) 

(sec-1 10-2)  (sec-110-2)  (sec-110-3) 

Rat  1  (Necrotic) 

1.45±0.013 

1.22±0.019 

1.86±0.017 

2.02±0.026 

2.74±0.041 

2.41±0.051 

4.05±0.059 

Rat  2  (Necrotic) 

3.48±0.048 

2.77±0.034 

4.28±0.048 

4.33±0.040 

2.98±0.048 

3.03±0.061 

4.76±0.062 

Rat  3  (Edematous) 

4.94±0.052 

5.16±0.067 

4.22±0.052 

4.13±0.067 

4.14±0.070 

4.27±0.078 

5.39±0.085 

Rat  4  (Edematous) 

5.25±0.053 

5.31  ±0.063 

5.07±0.068 

5.22±0.063 

4.43±0.075 

4.03±0.072 

3.85±0.056 

TABLE  IV 


Three-compartment  Model:  Estimated  pharmacokinetic  parameters  using  EKF  algorithm 


*.(3) 

roa 

(sec-110-2) 

kP 

(sec-1 10-2) 

1.(3) 

roc 

(sec-1 10-2) 

(sec-1 10-2) 

1.(3) 

Kout 

(sec-1 10-3) 

Rat  1  (Necrotic) 

1.93±0.061 

1.28±0.049 

1.82±0.032 

2.02±0.041 

3.89±0.052 

Rat  2  (Necrotic) 

4.41±0.074 

2.48±0.067 

4.87±0.066 

5.03±0.057 

5.45±0.071 

Rat  3  (Edematous) 

4.71±0.085 

3.88±0.077 

4.95±0.059 

4.68±0.050 

4.42±0.040 

Rat  4  (Edematous) 

5.29±0.091 

6.48±0.096 

4.48±0.062 

4.20±0.048 

5.01±0.055 

models  indicate  that  the  exchange  rates  between  the  capillary  and  the  adjacent  compartment,  k™,  k%, 
n  =  2,3,4,  are  significantly  different  for  the  necrotic  and  edematous  tissue.  We  observe  that  for  the 
four-  and  three-compartment  models,  the  estimated  exchange  rates  between  the  ISF  and  parenchymal 
cell  compartments,  A$,  n  =  3, 4,  are  comparable.  Similarly,  the  estimated  rate  of  drainage  out  of  the 
plasma,  k™ut,  n  =  2,3, 4,  are  consistent  for  all  models. 

Based  on  the  model  parameter  estimates,  we  computed  the  BIC  values  for  each  rat  data  to  reveal 
overfitting.  The  BIC  suggests  that  the  two-compartment  model  is  sufficient  for  all  four  measurement  sets. 
We  further  analyze  the  goodness-of-fit  of  the  compartmental  models  by  means  of  the  residual  analysis. 
The  details  and  results  explaining  why  two-compartment  model  is  a  better  statistical  fit  is  given  in  [5]. 
Figure  19  shows  the  measured  total  concentration  data  and  its  1-step  ahead  prediction  based  on  the 
two-compartment  model  for  each  rat  data.  Clearly,  there  is  a  good  agreement  between  the  actual  and  the 
estimated  measurements. 

As  we  concluded  that  the  two-compartment  model  provides  the  best  statistical  fit  for  the  rat  data,  we 
investigated  the  estimated  model  parameters  in  more  detail. 


ICG  concentration  (pm)  ICG  concentration  (pm) 
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TABLE  V 

Two-compartment  Model:Estimated  pharmacokinetic  parameters  and  volume  fractions  using  EKF 

ALGORITHM 


1.(2) 

A >a 

(sec-1 10-2) 

(sec-1 10-2) 

1.(2) 

" out 

(sec-1 10-3) 

(2) 

v\ 

(10-2) 

Vp2) 

(10-2) 

Rat  1  (Necrotic) 

2.47±0.043 

1.06±0.052 

4.61  ±0.073 

21.8±1.92 

1.41±0.053 

Rat  2  (Necrotic) 

3.54±0.082 

2.98±0.086 

4.83±0.092 

25.4±3.49 

2.42±0.088 

Rat  3  (Edematous) 

6.90±0.101 

4.93±0.072 

3.95±0.048 

30.4±2.81 

4.84±0.120 

Rat  4  (Edematous) 

8.40±0.114 

7.77±0.091 

4.02±0.068 

53.0±4.73 

7.03±0.321 

0  100  200  300  400  500  0  100  200  300  400  500 


Fig.  19.  ICG  concentration  measurement  data  and  1-step  prediction  of  the  measurements  for  four  different  rats. 
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In  order  to  characterize  the  difference  between  these  two  tumors,  we  estimated  the  pharmacokinetic 
parameters  ka,  kb  and  kout,  and  the  volume  fractions  vp  and  v j  for  each  rat  in  the  two-compartment 
model.  Table  V  tabulates  the  estimated  parameters.  The  rate  of  leakage  into  the  ISF  from  the  capillary, 
ka,  range  from  0.0247  to  0.0840  sec-1  and  the  rate  of  drainage  out  of  the  ISF  and  into  the  capillary,  kb, 
range  from  0.0106  to  0.0777  sec-1.  Note  that  the  permeability  rates  for  the  necrotic  cases  are  lower  than 
the  ones  observed  for  the  edematous  cases.  Additionally,  the  estimated  values  for  the  pharmacokinetic 
rates  are  much  higher  than  the  normal  tissue  values  due  to  the  increased  leakiness  of  the  blood  vessels 
around  the  tumor  region.  The  estimated  plasma  volume  fractions  agrees  with  the  values  reported  earlier 
[10].  These  results  confirm  that  vp  can  be  significantly  large  in  tumors  and  that  its  magnitude  varies  with 
respect  to  the  stage  of  the  tumor.  The  estimated  values  of  ISF  volume  fraction,  V{,  range  from  0.218  to 
0.53,  in  agreement  with  0.2  to  0.5  range  reported  earlier.  Note  that  these  results  are  valid  only  for  the 
ICG  pharmacokinetics  in  tumor  cells  R3230ac,  adenocarcinoma  and  may  not  be  generalized  for  other 
types  of  contrast  agents  or  tumor  types. 

Figure  20  shows  the  estimated  ICG  concentrations  in  the  plasma  and  the  ISF  compartments  for  the 
two-compartment  model  for  Rats  1  to  4.  The  peak  values  of  the  plasma  concentration,  Cp,  range  from 
2.72  M  to  4.28  ;iM.  The  absolute  value  of  the  concentrations  may  not  be  very  useful.  However, 
concentration  of  ICG  in  one  compartment  relative  to  the  concentration  in  another  compartment  may 
provide  useful  information.  We  consider  the  ratio  of  the  peak  concentrations  in  the  plasma  and  ISF  as 
a  potential  parameter  to  discriminate  different  tumors.  The  peak  Cp/Ci  ratio  for  Rats  1  to  4  is  0.551, 
0.593,  0.787,  1.151,  respectively.  This  ratio  is  higher  in  edematous  cases  consistent  with  the  fact  that 
ICG-albumin  leaks  more  into  the  ISF  in  advanced  tumors.  Additionally,  the  ICG  concentration  in  plasma 
decays  faster  than  the  ICG  concentration  in  ISF  due  to  ICG  elimination  through  the  liver  and  kidneys. 

2)  Breast  Data  Analysis:  In  order  to  further  investigate  the  pharmacokinetic  modeling  of  ICG  kinetics 
we  are  now  working  on  the  data  collected  from  three  different  patients  with  different  breast  tumor 
types.  First  case,  Case  1,  is  fibroadenoma,  which  corresponds  to  a  mass  estimated  to  be  1—2  cm  in 
diameter  and  located  around  6—7  o’clock  by  palpation  within  a  breast  of  9  cm  diameter.  Second  case. 
Case  2,  is  adenocarcinoma  corresponding  to  a  tumor  estimated  to  be  2—3  cm  in  diameter  and  located 
around  4—5  o’clock  by  palpation  within  a  breast  of  7.7  cm  diameter.  The  third  case,  Case  3,  is  invasive 
ductal  carcinoma,  which  corresponds  to  a  mass  estimated  to  be  3— 4  cm  in  diameter  and  located  around 
6  o’clock.  ICG  was  injected  intravenously  by  bolus  with  a  concentration  of  0.25  mg  per  kg  of  body 
weight.  Diagnostic  information  for  tumors  are  obtained  using  biopsy  results.  Since  biopsy  modifies  the 
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(a) 


(c) 


(b) 


Fig.  20.  ICG  concentrations  in  plasma,  Cp(t)  and  ISF,  Ci(t),  for  four  different  rats,  (a)  Rati,  (b)  Rat2,  (c)  Rat3,  (d)  Rat4. 

blood  volume  and  blood  flow  around  the  tumor  region,  measurements  were  made  before  the  biopsy.  Data 
acquisition  started  before  injection  of  ICG  and  continued  10  minutes. 

In  this  study,  we  make  use  of  the  data  collected  with  a  continuous  wave  (CW)  NIR  imaging  apparatus. 
The  apparatus  has  16  light  sources,  namely,  tungsten  bulbs  with  less  than  1  watt  of  output  energy.  They  are 
located  on  a  circular  holder  at  an  equal  distance  from  each  other  with  22.5  degree  apart.  Sixteen  detectors, 
namely,  silicon  photodiodes,  are  situated  in  the  same  plane.  The  breast  is  arranged  in  a  pendular  geometry 
with  the  source-detector  probes  gently  touching  its  surface.  Figure  21  illustrates  the  configuration  of  the 
apparatus  and  the  configuration  of  the  detectors  and  the  sources  in  a  circular  plane.  A  band  pass  filter  at 
805nm,  the  absorption  peak  of  ICG,  is  placed  in  front  of  the  sources  to  select  the  desired  wavelength. 
A  set  of  data  for  one  source  is  collected  every  500  ms.  The  total  time  for  a  whole  scan  of  the  breast 
including  16  sources  and  16  detectors  is  8.8  seconds.  The  detectors  use  the  same  positions  as  the  sources 
to  collect  the  light  originating  from  one  source  at  a  time.  Only  the  signals  from  the  farthest  1 1  detectors 
are  used  in  the  analysis.  For  example,  when  Source  1  is  on,  the  data  is  collected  using  Detectors  4  to 
14.  The  details  of  the  data  collection  procedure  and  the  apparatus  is  given  in  [23]. 
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Fig.  21.  The  cut  section  of  the  16  light  source-detector  device,  holding  a  human  breast  inside.  The  diameter  can  be  fitted 
easily.  The  16  light  source-detector  combinations  in  each  arm  are  located  equal  distance  (11  degree  apart),  but  when  the  device 
fits  the  breast,  only  the  diameter  chances. 

Using  the  CW  imager  described  above,  for  each  patient,  sufficient  number  of  source  detector  readings 
were  collected  from  different  angles.  Intes  et  al .  [24]  reconstructed  the  differential  absorption  images  using 
these  detector  source  readings  for  46  different  time  instants.  The  details  of  the  reconstruction  algorithm 
is  given  in  [24].  A  sample  set  of  reconstructed  differential  absorption  images  for  Case  1,  for  9  selected 
time  instants  are  given  in  Figure  22. 

In  Figure  22  each  one  of  9  reconstructions  consist  of  649  voxels  and  each  of  these  voxel  values 
correspond  to  a  differential  absorption  coefficient  value.  Using  the  the  linear  relationship  between  ICG 
concentration  and  absorption  coefficient,  we  get  ICG  concentration  maps  from  differential  absorption 
mappings  given  in  Figure  22.  A  sample  set  of  ICG  concentration  maps  for  the  selected  time  instants  is 
illustrated  in  Figure  23.  Here  the  concentration  maps  represents  bulk  ICG  concentration  in  the  tissue,  not 
specifically  in  the  plasma  or  the  ISF. 

ICG  concentration  curves  for  each  voxel  is  traced  using  all  reconstructions  at  different  time  intervals. 
For  Case  1,  we  have  46  reconstructions.  We  have  649  ICG  concentration  curves  representing  649  voxels 
and  each  of  these  curves  has  46  data  points  representing  46  different  time  points.  An  example  of  the 
ICG  concentration  curve  for  the  voxel  where  maximum  uptake  has  occurred  is  given  in  Figure  24. 

Using  two  compartmental  model ,  we  estimated  pharmacokinetic  parameters  for  each  voxel  using  EKF 
algorithm.  We  then  constructed  2-D  permeability  rate  mappings  using  values  of  these  parameters.  2-D 
maps  for  ka,  k ^  for  two-compartmental  model  are  presented  in  Figures  25(a),  25(b). 
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Fig.  22.  Differential  absorption  reconstruction  maps  for  the  9  different  time  instances. 
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Fig.  24.  ICG  concentration  curve  for  a  specific  voxel  obtained  using  46  different  reconstructions  at  different  time  instants. 
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(a)  2-D  spatial  map  of  the  pharmacokinetic 
rates  representing  the  drainage  outside  of 
the  plasma 
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(b)  2-D  spatial  map  of  the  pharmacoki¬ 
netic  rates  representing  the  leakage  into  the 
plasma 


Fig.  25.  2-D  spatial  mapping  of  pharmacokinetic  rates  representing  the  drainage  outside  of  the  plasma  (a)  and  into  the  plasma 

(b) 
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We  then  analyzed  if  the  kinetic  rates  are  statistically  different  for  inside  and  outside  of  the  tumor 
region.  The  ka  and  kf,  values  from  inside  and  outside  the  tumor  region  are  statistically  different  with 
a  p-value  of  0.0001.  Using  these  results  one  may  conclude  that  spatially  resolved  pharmacokinetic  rate 
images  can  be  used  for  breast  cancer  screening  and  diagnosis. 

IV.  Key  Research  Accomplishments 

The  SOW  tasks  defined  for  the  first  12  months  of  the  project  have  been  completed.  In  particular, 

1)  An  MR  image  analysis  toolbox  was  developed.  The  toolbox  performs  image  segmentation,  extrac¬ 
tion  of  anatomical  prior  information  and  co-registration  of  anatomical  and  optical  images.  (Aim  1 , 
Task  1). 

2)  Two-level  domain  decomposition  methods  has  been  developed  for  diffuse  optical  tomography.  The 
convergence  of  these  methods  are  shown.  The  methods  developed  give  smaller  mean  square  error 
with  less  computational  cost  than  the  methods  without  domain  decomposition  techniques.  The 
reconstruction  algorithms  can  be  made  more  efficient  with  parallel  computing.  If  the  tumor  region 
of  the  optical  coefficients  is  known  by  a  priori  anatomical  information,  the  reconstruction  algorithms 
become  much  more  efficient  and  accurate  by  updating  around  the  tumor  region  only.  (Aim  1,  Task 
2) 

3)  The  hierarchical  Bayesian  approach  has  been  successfully  applied  for  the  recovery  of  absorption 
coefficients  [4]  and  for  the  direct  recovery  of  tissue  chromophores  [2]  and  has  been  proven  to  be 
effective  in  addressing  the  correlation  between  anatomical  and  functional  images  (Aim  1,  Task  3 
and  4) 

4)  Three  different  compartmental  models  have  been  developed  for  the  flow  kinetics  of  the  optical 
contrast  agent  ICG.  A  robust  and  efficient  method  that  can  simultaneously  estimate  permeability 
rates  and  contrast  agent  kinetics  in  various  compartments  was  developed.  (Aim  2,  Task  1  and  2) 

5)  The  computational  models  and  methods  developed  were  applied  to  Fisher  rat  data  carrying  adeno¬ 
carcinoma  and  shown  that  two  compartment  model  is  sufficient  to  characterize  the  kinetics  of  the 
ICG  and  to  differentiate  different  tumor  types.  (Aim  2,  Task  2  and  3) 

6)  A  method  to  construct  spatially  varying  pharmacokinetic  rates  is  developed  and  is  currently  being 
applied  to  human  breast  data.  (Aim  2,  Task  3) 
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V.  Reportable  Outcomes 

The  hierarchical  Bayesian  approach  has  been  successfully  applied  for  the  recovery  of  absorption 
coefficients  [4]  and  for  the  direct  recovery  of  tissue  chromophores  [2]. 

Two-level  multiplicative  overlapping  domain  decomposition  method  (TMODDM)  is  applied  to  the 
linearized  DOT  [52],  With  TMODDM  as  a  forward  solver,  two-level  multiplicative  space  decomposition 
method  is  adopted  to  the  nonlinear  DOT  in  [53],  [54],  [55].  The  convergence  and  computational  costs  of 
the  proposed  algorithms  are  shown  in  [55].  The  reconstruction  of  diffusion  coefficient  is  shown  in  [53] 
and  the  reconstruction  of  absorption  coefficients  is  shown  in  [54],  [55]. 

Compartmental  modeling  of  ICG  pharmacokinetics  using  a  two-compartment  model  for  Fisher  rat 
data  has  presented  in  [6],  Compartmental  modeling  of  ICG  pharmacokinetics  using  four,  three,  and  two- 
compartment  models  and  the  best  model  order  selection  criteria  has  presented  in  [5]  as  a  journal  paper. 
A  summary  of  this  work  has  also  presented  in  [6]  and  in  [8].  Spatial  mapping  of  pharmacokinetic  rates 
using  breast  tumors  have  been  presented  in  [9]  Era  of  Hope  meeting  for  the  Department  of  Defense 
(DOD)  Breast  Cancer  Research  Program  in  June,  2005. 

Complete  list  of  reportable  outcomes  is  given  below: 

1)  Intes  X,  Maloux  C,  Guven  M,  Yazici  B  and  Chance  B  2004  Diffuse  optical  tomography  with 
physiological  and  spatial  a  priori  constraints  Phys.  Med.  Biol.  49  N 155-63 

2)  Guven  M,  Yazici  B,  Intes  X  and  Chance  B  2005  Diffuse  optical  tomography  with  a  priori  infor¬ 
mation  Phys.  Med.  Biol.  50  2837-58 

3)  M.  Guven,  B.  Yazici,  X.  Intes,  B.  Chance,  “3D  Diffuse  Optical  Tomography  with  a  priori  Anatom¬ 
ical  Information”  SPIE-OSA  Joint  Proc.  SPIE  Vol.  5138,  p.  268-280. 

4)  M.  Guven,  B.  Yazici,  X.  Intes,  B.  Chance,  2003,  “An  Adaptive  V-grid  algorithm  for  Diffuse  Optical 
Tomography,”  International  Conference  on  Image  Processing,  2,  14-7. 

5)  M.  Guven,  B.  Yazici,  X.  Intes,  B.  Chance,  “Recursive  least  squares  algorithm  for  optical  diffusion 
tomography,”  2002,  Bioengineering  Conference  Proceedings  of  the  IEEE  28th  Annual  Northeast, 
273-4. 

6)  B.  Alacam,  B.  Yazici,  X.  Intes,  B.  Chance  ”  Extended  Kalman  Filtering  for  the  Modeling  and 
Analysis  of  ICG  Pharmacokinetics  in  Cancerous  Tumors  using  NIR  Optical  Methods”,  Transactions 
in  IEEE  Biomedical  Engineering ,  In  Review. 

7)  B.  Alacam,  B.  Yazici,  X.  Intes,  B.  Chance  ’’Extended  Kalman  Filtering  Framework  for  the  Modeling 
and  Analysis  of  ICG  Pharmacokinetics,”  Proc.  Of  2005  SPIE  Photonic  West,  San  Jose,  California 
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USA,  22  -  27  January  2005,  vol.  5693.  pp.  17-27. 

8)  B.  Alacam,  B.  Yazici,  X.  Intes,  B.  Chance  ’’Analysis  of  ICG  Pharmacokinetics  in  Cancerous  Tumors 
using  NIR  Optical  Methods,”  Proc.  of  EMBS  —  27th  Anniversary  Conference,  Shanghai,  China, 
September  2005. 

9)  B.  Yazici,  K.  Kwon,  M.  Guven,  B.  Alacam,  ”  Concurrent  MR-NIR  Imaging  for  Breast  Cancer 
Diagnosis”,  The  Breast  Cancer  Imaging  Network  NTROI  Annual  Network  Retreat  and  Workshop, 
Newport  Beach,  June  2-4,  2005. 

10)  B.  Yazici,  B.  Alacam,  X.  Intes,  N.  Shoko,  B.  Chance,  ”  Compartmental  Modeling  of  ICG  Pharma¬ 
cokinetics  for  Breast  Cancer  Tumors”  Era  of  Hope  meeting  for  the  Department  of  Defense  (DOD) 
Breast  Cancer  Research  Program,  Philadelphia,  June  8-11,  2005. 

11)  Il-Y.  Son,  M.  Guven,  B.  Yazici,  and  X.  Intes,  2004  A  2-level  domain  decomposition  algorithm  for 
inverse  diffuse  optical  tomography,  International  Conference  on  Image  Processing,  Proceedings  of 
ICIP,  5,  3315  -  3318. 

12)  K.  Kwon,  I.-Y.  Son,  and  B.  Yazici,  2005  Domain  decomposition  method  for  diffuse  optical  tomog¬ 
raphy,  Computational  Imaging  III ,  edited  by  Charles  A.  Bouman  and  Eric  L.  Miller,  Proceedings 
of  SPIE  ,  volume  5674,  64-75. 

13)  K.  Kwon,  I.-Y.  Son,  and  B.  Yazici,  2005  Two-level  overlapping  domain  decomposition  algorithm 
for  a  nonlinear  inverse  DOT  problem,  Optical  Tomography  and  Spectroscopy  of  Tissue  VII ,  edited 
by  B.  Chance,  R.  R.  Alfano,  B.  J.  Tromberg,  M.  Tamura,  and  E.  M.  Sevick-Muraca,  Proceedings 
of  SPIE  ,  volume  5693,  459^168. 

14)  K.  Kwon  and  B.  Yazici,  2005  Two-level  domain  decomposition  methods  for  diffuse  optical  tomog¬ 
raphy,  in  preparation. 


VI.  Conclusions 

In  the  last  12  months,  the  research  project  has  progressed  as  planned.  Our  research  demonstrated 
that  the  computational  requirements  and  accuracy  of  diffuse  optical  tomography  can  be  significantly 
improved  by  the  incorporation  of  a  priori  anatomical  information  provided  by  MRI.  We  developed  new 
compartmental  models  and  estimation  methods  to  analyze  the  flow  kinetics  of  the  optical  contrast  agent 
ICG  and  applied  our  approach  to  Fisher  rat  data  containing  adenocarcinoma.  Our  results  indicate  that  the 
flow  kinetic  parameters  can  be  used  to  distinguish  different  tumor  types.  Currently,  we  are  developing 
methods  to  build  spatially  resolved  flow  kinetic  images  using  high  temporal  resolution  NIR  spectroscopy 
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data.  We  are  in  the  process  of  applying  these  methods  to  human  breast  data.  These  results  were  presented 
in  2005  Era  of  Hope  Meeting,  Philadelphia,  PA.  Presently,  we  expect  the  next  12  months  of  research  to 
progress  in  line  with  the  tasks  outlined  in  the  SOW. 
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Abstract.  Diffuse  optical  tomography  (DOT)  in  the  near  infrared  involves  the 
reconstruction  of  spatially  varying  optical  properties  of  turbid  medium  from  boundary 
measurements  based  on  a  forward  model  of  photon  propagation.  Due  to  the  highly 
non-linear  nature  of  DOT,  high  quality  image  reconstruction  is  a  computationally 
demanding  problem  that  requires  repeated  solutions  of  both  the  forward  and  the 
inverse  problems.  Therefore,  it  is  highly  desirable  to  develop  methods  and  algorithms 
that  are  computationally  efficient.  In  this  paper,  domain  decomposition  methods 
are  introduced  to  address  the  computational  complexity  of  the  DOT  problem.  A 
two-level  multiplicative  overlapping  domain  decomposition  method  for  the  forward 
problem  and  a  two-level  multiplicative  space  decomposition  method  for  the  inverse 
problem  are  developed.  We  show  the  local  convergence  of  the  inverse  solver  and  derive 
the  computational  complexity  of  each  method.  The  performance  of  the  proposed 
methods  in  numerical  simulations  are  demonstrated.  The  convergence  analysis  and 
computational  cost  analysis  of  both  forward  and  inverse  solvers  are  provided. 
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1.  Introduction 

Diffuse  optical  image  reconstruction  based  on  the  diffusion  equation  is  a  highly  non¬ 
linear  ill-posed  problem  that  calls  for  the  use  of  non-linear  minimization  methods  with 
regularization  to  stabilize  the  solution  [1]. 

DOT  reconstruction  is  posed  as  an  optimization  problem  involving  two  coupled 
steps,  namely  forward  and  inverse  problems.  Each  consists  of  an  iterative  solver  whose 
solutions  are  used  as  an  input  to  the  other  solver.  More  precisely,  the  forward  solver 
computes  the  photon  density  and  its  Jacobian  with  respect  to  the  optical  coefficients, 
and  the  inverse  solver  updates  the  optical  coefficients  based  on  the  output  of  the  forward 
step.  The  updated  coefficients  are  then  used  in  the  forward  solver  to  recompute  the 
photon  density  and  its  Jacobian.  As  a  result,  the  computational  complexity  of  the 
DOT  reconstruction  using  these  approaches,  quickly  grows  with  the  number  of  pixels 
and  dimensions.  Thus,  real  time  computation  of  DOT  requires  numerical  techniques  to 
reduce  the  complexity  of  the  problem. 

In  this  paper,  we  investigate  domain  decomposition  methods(DDM)  to  address  the 
computational  requirements  of  the  DOT  image  reconstruction.  Domain  decomposition 
methods  have  been  developed  in  the  last  two  decades  in  the  area  of  numerical  solutions 
of  partial  differential  equations  motivated  by  the  need  for  fast  and  efficient  algorithms  for 
solving  large-scale,  three-dimensional  problems.  A  major  advantage  of  the  DDM  is  that 
they  allow  for  distributed  parallel  numerical  solvers  on  smaller  subdomains,  making  the 
computation  extremely  efficient.  We  applied  DDM  to  the  forward  problem  formulated 
by  the  solution  of  an  elliptic  partial  equation  given  optical  coefficients  and  space 
decomposition  method  to  the  inverse  problem  formulated  as  a  minimization  problem 
given  boundary  measurements.  Badea  et  al.  introduced  the  term  space  decomposition 
to  describe  the  application  of  DDM  approach  to  the  constrained  convex  minimization 
problems  arising  from  variational  inequalities  and  provided  conditions  under  which  DDM 
based  approaches  lead  to  the  convergence  [4] -  [6] ,  [23]  -  [25] .  DDM  begins  by  partitioning 
the  domain  into  two  or  more  subdomains  as  in  Figure  1.  The  inverse  problem  is  then 
divided  into  subproblems  on  each  of  the  subdomains. 

Among  many  kinds  of  DDM  techniques  we  have  considered  the  multiplicative 
overlapping  domain  decomposition  method  (MODDM)  for  the  forward  problem  and 
multiplicative  space  decomposition  (MSDM)  for  the  inverse  problem.  These  methods  are 
called  one-level  methods.  The  proposed  methods  in  this  paper  are  two-level  MODDM 
(TMODDM)  and  two-level  MSDM  (TMSDM)  for  forward  and  inverse  problems, 
respectively. 

Since  the  photon  density  function  is  rapidly  varying  around  the  point  source,  the 
contribution  from  the  sources  may  be  neglected  in  the  subdomains  that  are  far  apart 
from  the  sources.  Thus  the  one-level  MODDM  is  not  sufficient  to  approximate  the 
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photon  density  function  for  the  point  source,  since  the  subdomain  correction  step  is 
confined  to  each  subdomain.  Therefore  we  propose  an  algorithm  based  on  TMODDM 
for  the  forward  solver  [21].  The  coarse  level  correction  step  applied  to  the  entire  domain 
decreases  low  frequency  errors  that  are  not  handled  well  at  the  subdomain  correction 
step.  Therefore  TMODDM  decreases  the  overall  error  more  efficiently  than  MODDM. 

The  use  of  two-level  method  is  extended  to  the  inverse  problem.  The  convergence  of 
MSDM  is  shown  in  [5]  under  the  some  conditions  for  the  objective  function  F.  Among 
the  conditions,  the  strong  convexity  condition  for  F,  that  is  the  positive-definiteness  of 
the  Hessian  of  F,  can  be  shown  only  locally.  Thus  to  prove  the  convergence  of  MSDM, 
we  require  that  the  optical  coefficients  for  each  iteration  step  and  at  each  subdomain 
must  be  sufficiently  close  to  the  target  optical  coefficient.  In  other  words,  a  good  initial 
guess  for  the  target  optical  coefficient  is  required.  Therefore,  we  have  chosen  upsampled 
coarse  level  image  as  an  initial  guess  for  Algorithm  MSDM  assuming  it  is  sufficiently 
close  at  least  at  low  frequency. 

The  paper  is  organized  as  follows:  In  the  next  section,  the  forward  solver,  based  on 
the  photon  diffusion  equation  model  and  minimization  formulation  of  the  inverse  DOT, 
is  described.  Algorithms  applying  domain  decomposition  techniques  to  the  forward  and 
inverse  steps  of  the  minimization  formulation  are  presented  in  Section  3.  Section  4 
presents  numerical  simulations  and  Section  5  summarizes  our  results  and  conclusions. 
The  paper  includes  an  appendix  for  the  proof  of  the  local  convergence  of  TMSDM 
algorithm. 

2.  Diffuse  Optical  Tomography 

2.1.  Photon  Diffusion  Equation  in  Frequency  Domain 

Propagation  of  light  in  biological  tissues  is  well-modelled  by  the  photon  diffusion 
equation  with  the  Robin  boundary  condition.  In  frequency  domain  photon  diffusion 
equation  is  given  as  follows: 

qfj) 

-V  •  (kV$)  +  (pa  +  —  )$  =  q  in  fi 

c 

$  +  2<m—  =  0  on  dfl  (2.1) 

ou 

where  f2  is  a  Lipschitz  domain  in  R",  n  =  2,3,  dfl  is  its  boundary,  u  is  the  unit  outward 
normal  vector  on  the  boundary,  <f>  is  the  photon  density,  q  is  a  source  term,  and  pa,  pfs, 
and  k  =  3^- y  are  the  absorption,  reduced  scattering,  and  diffusion  coefficients, 
respectively.  The  constant  R,  is  determined  by  the  refraction  on  the  boundary. 

The  unique  identification  of  the  optical  coefficients  pa  and  k  in  (2.1)  when  Dirichlet- 
to-Neumann  map  is  given  (or  infinite  sources  and  infinite  detectors  are  given),  can  be 
easily  shown  by  using  the  uniqueness  results  for  isotropic  case  [22] .  For  the  uniqueness 
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of  the  optical  coefficients  when  k  has  anisotropic  anomalous  region  contained  in  a  known 
background,  see  [10,  11,  12]. 

On  the  boundary,  we  measure  the  following  Rytov  data  [18]: 


r 


=  ^  (2-2a) 

=  los(^)'  <2-2b> 


[2,  3]  show  numerical  examples  where  the  Rytov  approximation  provides  much  better 
images  than  the  Born  approximation. 

In  this  paper,  we  assume  that  p's  is  known  and  concentrate  only  on  the 
reconstruction  of  the  absorption  coefficient  pa.  Thus,  the  Jacobian  of  T  with  respect  to 
pa  for  given  p's  as  follows: 

■J— -(r)  =  [  [-3K2VG(r,r,)vW)  +  G(r,r,)¥{?)}dr,,  (2.3) 

dpa  2 ad>  Jn 

where  $  is  the  complex  conjugate  of  d>. 


2.2.  Discretization  and  Finite  Element  Method 

Suppose  Ns  sources  r^,  j  =  1,  •  •  • ,  Ns  and  Nd  detectors  ri+Ns,  i  —  1,  •  •  • ,  Nd  on  boundary 
of  ft  is  given.  Let  d>j  be  the  solution  of  (2.1)  for  point  source  qj(r)  =  6(r  —  rf),j  = 
1,  •  •  • ,  Ns.  Define 

r<j  =  log(^(r<+Jv.)),  (2.4) 

for  i- th  detector  and  j-th  source. 

Consider  the  finite  element  space  with  Nn  Uk,k  =  1,  •  •  • ,  Nn  bases.  For  piecewise 
bilinear  element,  Nn  is  the  same  as  the  number  of  nodes.  Let  the  number  of  elements 
be  Ne  and  denote  the  elements  by  Tm,  m  =  1,  •  •  • ,  Ne. 

The  finite  element  formulation  for  (2.1)  in  this  finite  element  space  for  the  j-th. 
source  qj  (j  —  1,  •  •  • ,  Ns)  is  as  follows: 

K  +  C+^-A 
2  a 

where  K,  C,  A  are  Nn  x  Nn  matrices  and  IP  is  a  Nn  x  1  column  vector  given  by 

Kki  =  /  KVukf7ui 

Ju 


(2.5) 
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V  = 


1 

0 


if  l  =  j 

Otherwise, 


for  k,  l  =  1,  •  ■  • ,  Nn. 

Assume  that  L  and  U  are  constants  such  that  0  <  L  <U.  Define  a  function  space 
for  pa  as 

V  =  {Pa  e  L2(ty\L  <  \pa(r)\ren  <  U}, 

Vxe  —  {pa  G  V\ pa  is  constant  at  each  Tm,m=  1,  ■  •  • ,  Ne}. 

Define  norms  in  V  and  VNe  as  follows: 


x 


=  for  x  e  V> 


\x\\vNe  - 


Ne 

E 


x(my 


\  tx  lT” 


for  x  G  VNe, 


where  x(m)  is  the  constant  at  the  element  Tm  and  \Tm\  is  the  area  of  Tm,  m  =  1,  •  •  • ,  Ne. 

Assume  that  k,  pa  G  V^e.  Let  ^j(k)  be  the  value  of  at  fc-th  node  point  and 
/c(m )  be  the  value  of  k  at  m-th  element  Tm.  Prom  the  solution  of  (2.5),  we  obtain  the 
boundary  measurement  data  By  discretizing  (2.3),  we  get  the  value  of  the  Jacobian 


of  at  the  m-th  element: 
dTi 


Nn 


g^(m)  =  "  E  l)  +  Wm(k,  !)]*,(«) 


Vm 

W, 


(k,l)  =  / 
J  Tn 

n(k,l)  =  j 

J  Tn 


WukVui, 


ukut. 


(2.6) 


2.3.  DOT  as  a  nonlinear  ill-posed  optimization  problem  and  the  trust  region  method 

Given  (2.1)  and  (2.6),  DOT  is  formulated  as  the  following  nonlinear  minimization 
problem  to  determine  the  absorption  coefficient  pa  in  Vnb  '■ 

min^a€VNeF(pa),  (2.7) 

1  Ns  Nd  Ne 

fm  =  2  E  E  w  -  m‘F + f  E 

j— 1  i= 1  m=  1 

where  Mij  are  measured  data  at  the  i-th  detector  and  the  j- th  source,  and  a  is  a 
regularization  parameter. 

Although  the  uniqueness  of  pa  (  and  p's )  is  known  for  infinite  sources  and  detectors, 
the  unique  solvability  of  (2.7)  for  finite  sources  and  detectors  is  not  known  leading  to  the 
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ill-posedness  of  the  DOT  problem.  To  address  this  ill-posedness,  nonlinear  optimization 
such  as  Newton  method  is  employed. 

The  optimization  problem  in  (2.7)  is  composed  of  two  steps;  the  step  to  determine 
the  minimizing  direction  at  the  current  coefficients  pa  and  the  step  to  perform  line  search 
on  those  minimizing  directions.  In  the  Newtonian  method,  the  minimizing  direction  is 
—{F")~lF' ,  where  F'  and  F"  are  the  Jacobian  and  the  Hessian  of  F  with  respect  to  pa. 
In  the  classical  Newton  method,  searching  direction  is  found  by  computing  F"  directly. 
However,  this  method  requires  extensive  computation  and  is  unstable.  To  decrease 
complexity  of  the  computation,  in  quasi-Newton  methods,  F"~l  is  approximated  by 
various  computationally  efficient  methods  such  as  the  Davidon-Fletcher-Powell  method, 
the  Broyden-Fletcher-Goldfarb-Shanno  method,  and  the  conjugate  gradient  method  [7]. 
To  stabilize  the  computation,  F"  +  XI  is  used  in  place  of  F"  in  the  Levenberg-Marquadt 
method  [16,  17]. 

In  this  paper,  we  adopted  the  trust  region  method,  a  variant  of  Newton  approach, 
as  an  optimizer.  The  trust  region  method  solves  the  following  quadratic  minimization 
problem: 

min  S^8pa  G  W,  \\D6 fia\\  <  T  ^(6pa)tF"(pa)5pa  +  (8paYF'(pa)  j  (2.8) 

where  D  is  a  scaling  matrix,  T  is  a  trust  region  parameter,  and  W  is  a  subspace  of  V. 
The  scaling  matrix  D  is  used  to  handle  constraints  for  the  minimization.  To  stabilize  the 
minimization,  we  control  the  trust  region  parameter  T,  which  is  similar  to  control  the 
A  parameter  in  the  Levenberg-Marquadt  method  [7].  To  avoid  extensive  computation, 
the  subspace  W  is  chosen  as  the  two-dimensional  subspace  composed  of  the  gradient 
direction  and  the  approximate  Newton  direction  [9].  This  makes  the  trust  region  method 
suitable  for  large-scale  constrained  optimization  problems  like  the  one  state  in  (2.7). 

Taking  vanishing  gradient  point  of  the  quadratic  form  ^{8pa)tF"{pa)8pa  + 
(SpaYF^Spa)  in  (2.8),  we  get 

( Jt{8pa)J(5pa )  +  aINe)  8 pa  =  -(«/(/*„)*&(/«„)  +  ocpa).  (2.9) 

Thus,  if  assume  that  8pa  is  sufficiently  small,  the  trust  region  method  can  be  used  to 
solve  (2.8)  at  each  iteration. 

3.  Two-level  Domain  decomposition  methods  for  the  diffuse  optical 
tomography 

In  this  section,  we  describe  the  two-level  domain  decomposition  methods  considered  in 
this  paper,  as  applied  to  the  forward  and  inverse  problems.  These  two  methods  are 
composed  of  the  coarse  level  and  the  subdomain  correction  steps.  We  will  describe  our 
notation  and  approach  for  the  two-dimensional  optical  domain  0  =  [a,  b\  x  [c,  d]  C  R2 
and  bilinear  finite  element.  Its  extension  to  three-dimensional  domain  is  straightforward. 
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Let  0NxtNy  be  0  uniformly  divided  by  Nx  times  in  the  z-axis  direction  and  Ny  times 
in  the  y- axis  direction.  Thus  0NxtNy  has  Nn  =  (Nx+ 1) x ( Ny+1 )  nodes  and  Ne  —  NxxNy 
elements.  We  call  0Nx/^Ny/2  the  coarse  level  of  0NxtNy  assuming  Nx  and  Ny  are  even. 
Let  ClNXtNy  be  decomposed  into  disjoint  union  of  d  subdomains  Op,p  =  1  such 

that 

d 

0 Nx,Ny  =  (J  Op.  (3.10) 

P=1 

(3.10)  describes  the  nonoverlapping  domain  decomposition.  For  the  overlapping  domain 
decomposition,  we  define  Of,p  =  1 ,  •  •  • ,  d.  an  extension  of  Op,  recursively  for  all 
nonnegative  integers  w  as  follows:  =  Op  and  Of  is  the  union  of  Of-1  and  its  adjacent 

elements  in  contact  with  the  boundary  of  Of*1,  where  w  will  be  called  the  overlapping 
width.  Thus  the  overlapping  domain  decomposition  is  given  by 

d 

=  (J  nr  ■  (3-11) 

P=1 


3.1.  Two-level  Multiplicative  Overlapping  Domain  Decomposition  Method  for  the 
forward  problem 

In  the  forward  solution  of  the  DOT  problem,  one  has  to  compute  the  photon  densities 
satisfying  (2.1)  and  its  Jacobian.  Commonly  used  methods  are  the  finite  element  method 
in  (2.5)  and  the  adjoint  method  (2.6).  We  apply  TMODDM  to  (2.5)  in  order  to  reduce 
the  computation  complexity  of  the  problem.  TMODDM  composed  of  two  steps:  the 
coarse  level  correction  and  the  subdomain  correction  steps. 

In  the  subdomain  correction  step,  we  use  MODDM.  Let  be  the  n-th  step  solution 
of  TMODDM  and  <J>n+d+i  be  the  updates  via  coarse  level  correction  to  be  presented 
below.  Then  the  updates  <hn+ci+1  at  p-th  subdomain  Of  is  obtained  by 


,n+£±i  f  v  in  Of 

$n+d+i  =  J  p  p 

1  in  0  \TT£ 

(3.12) 

where  v  is  the  solution  of  the  following  partial  differential  equation: 

-V  •  (kVv)  +  (pa  +  —  )v  ==  q 

c 

ino; 

(3.13a) 

n  dv 

v  +  2  an—  —  0 
ov 

on  dOf  fl  dO, 

(3.13b) 

V  =  $n+dfr 

on  dOf  \  dO. 

(3.13c) 

We  call  the  coarse  level  correction  step  combined  with  MODDM  TMODDM.  In  the 
coarse  level  correction  step,  an  important  point  is  to  define  the  interpolation  operator 
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p  ■  ^Nx/2,Ny/2  -*  ^Nx,Ny  and  restriction  operator  R  :  0Nx<Ny  —>  0Nx/2.Ny/2  between  the 
functions  defined  at  the  node  points  in  £Inx/2,nv/2  and  0^XtNy. 

We  use  bilinear  interpolation  operator  P  and  full  weighting  restriction  operator  R, 
having  the  following  symbols: 

l  1/4  1/2  1/4  \ 

P  =  1/2  1  1/2  and  R  =  -PT, 

\  1/4  1/2  1/4  ) 

where  PT  is  the  transpose  of  P[8,  24].  At  the  boundary,  P  and  R  has  the  symbol 
(1/2  1  1/2)  boundary  is  a  one-dimensional  line. 

Let  the  maximum  sweep  of  all  subdomain  corrections  followed  by  the  coarse  grid 
correction  be  MF.  A  pseudo-code  summarizing  the  Algorithm  TMODDM  to  solve  (2.5) 
at  D,Nx,Ny  is  given  as  follows: 

Algorithm  1  TMODDM 
for  j  =  1  •  •  •  Ns  +  Nd  do 
Initialize 

for  n  =  1  ■  •  •  Mf  do 

rj  •<=  R((K  +  C  +  ^A)4>”-1  -  V),  rj  G  0Nx/2,Ny/2  {Restrict  the  residual.} 
es  4=  (K  +  C  +  ~A)“V,  ej,rj  G  DNx/2)Ny/2  {Coarse  grid  correction.} 

4>"-1+ra  <=  +  P(ej),  P(ej)  €  nN„Ny 

for  p  =  1,  •  •  • ,  d  do 

„  1 1  p+i 

Update  d+1  at  Of  by  (3.12)  {Subdomain  correction.} 

end  for 
end  for 
end  for 

Compute  Fij  and  ^/i(m), m  =  1  ,Ne  using  (2.4)  and  (2.6).  {Post-processing} 


3.2.  Two-level  multiplicative  space  decomposition  method  for  the  inverse  problem 

In  the  forward  problem,  we  obtain  the  value  of  photon  density  at  each  nodes.  Whereas 
in  the  inverse  problem,  we  are  interested  in  the  value  of  the  optical  coefficients  at  each 
element.  It  is  possible  to  define  different  domain  decomposition  schemes  for  the  nodes 
and  elements  in  the  forward  and  inverse  problems.  In  this  paper,  however,  we  use  the 
same  domain  decomposition  scheme  for  both  problems.  For  the  domain  decomposition 
(3.10)  and  (3.11),  let  the  restrictions  of  VJve  to  the  overlapping  subdomain  Of  and 
nonoverlapping  subdomain  Op  be  Vff  and  Wffe  for  p  =  1,  •  •  • ,  d,  respectively,  such  that 

Vjfe  =  {xe  VNe\x  =  0onQ\ttf},  (3.14a) 
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Wff  =  {x  G  Wnb \x  =  0  on  9  \  !+}.  (3.14b) 


TMSDM  composed  of  two  steps:  The  coarse  level  correction  and  the  subdomain 
correction  steps. 

The  subdomain  correction  step  uses  the  algorithm  MSDM  which  is  given  below:  Let 
pQa  be  the  initial  guess  for  optical  coefficient  on  fi nx,nv-  Define  pa+p'd  as  the  n-th  update 

of  the  optical  coefficient  in  the  p-th  subdomain.  Then  given  pa  d  6  .  p  =  1,  ■  ■  • ,  d, 

compute  5 pa  d  such  that 

5pa+d=argmin  n+§  n+E-i  n+*  F(pa+Ed~  +  5pa+d),  (3.15) 

Sva  deVpe,fla  d  +SfJ-a  d£VNe 

and  update 


n+ 


£ 

d 


(3.16) 


The  convergence  of  MSDM  is  proved  under  some  conditions  on  the  objective  function  F 
[5] .  We  have  proved  the  convergence  of  MSDM  in  Theorem  1  when  approximating  the 
pa  in  each  iteration  is  sufficiently  close  to  the  true  pa.  As  a  result,  the  convergence  of  the 
MSDM  is  assured  when  the  initial  guess  is  sufficiently  close  to  the  true  pa.  Motivated 
by  the  Algorithm  TMODDM,  we  applied  a  two-level  approach  to  the  MSDM  by  using 
the  coarse  level  solution  as  an  initial  guess.  For  the  present  work,  we  assume  that  the 
coarse  level  solution  is  sufficiently  close  to  the  true  image  of  the  absorption  coefficient. 

To  use  the  coarse  level  solution  in  the  finer  level  as  an  initial  guess,  we  adopted 
an  upsampling  operator  U  from  the  coarse  level  to  the  fine  level.  Let  y  be  a  function 
defined  on  Nx/2  x  Ny/2  elements  on  the  coarse  level,  then  U(y)  is  a  function  defined  on 
each  Nx  x  Ny  elements  on  the  fine  level  such  that 

U(y){2mx  -  1,2 my  -  1)  =  y(mx,my),  U(y)(2mx  -  l,2my)  =  y(mx,my), 
U(y)(2mx,2my  -  1)  =  y(mx,my),  U(y)(2mx,2my)  =  y(mx,my), 

for  1  <  mx  <  Nx/2, 1  <  my  <  Ny/2. 

Let  Me  and  Ms  denote  the  number  of  sweeps  for  the  coarse  and  subdomain 
corrections.  A  pseudo-code  summarizing  the  Algorithm  TMSDM  is  given  as  follows: 

Assume  that  we  know  the  location  of  the  anomalous  region  is  contained  in  some 
small  region  by  a  priori  information  provided  by  a  secondary  imaging  modality  such  as 
MRI  or  X-rays  or  by  a  posteriori  information  after  the  coarse  level  correction  step.  Then 
instead  of  performing  the  subdomain  correction  step  on  all  of  the  subdomains,  we  can 
update  pa  on  a  region  of  anomaly.  The  convergence  of  the  algorithm  when  implemented 
only  for  a  region  of  interest  does  not  change  if  the  initial  guess  of  the  optical  coefficient 
is  sufficiently  close  to  the  true  image  of  the  coefficient. 
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Algorithm  2  TMSDM 

Tc  =  log(R(exp(T)))  {Restrict  boundary  measurements.} 
for  i  =  1,  •  •  • ,  Me  do 

Solve  (2.8)  on  the  coarse  level  from  the  boundary  measurement  data  T?  •  to 
approximate  pca  {Coarse  level  correction.} 
end  for 
Ta  <=  U(pCa) 
for  p  =  1,  •  •  • ,  d  do 
for  k  =  1,  •  •  • ,  Ms  do 

Solve  (3.15)  from  the  measurement  data  Ty  and  elaborate  pa  on  0,'f  by  (3.16) 
using  Algorithm  1  as  a  forward  solver  {Subdomain  correction.} 

end  for 
end  for 

3.3.  Convergence  of  TMODDM  and  TMSDM 

In  this  subsection,  we  will  discuss  the  convergence  behavior  of  the  proposed  algorithms. 
Assume  that  the  mesh  size  of  the  finite  element  formulation  is  0(h)  and  that  the 
subdomains  are  of  diameter  0(H)  and  the  overlap  region  is  width  0(6H).  Then  the 
following  convergence  behavior  is  known  for  the  Algorithm  TMODDM  [20,  19]. 

•  Convergence  is  poor  if  8  =  0  but  improves  rapidly  as  5  increase.  (3.17a) 

•  If  5  is  fixed,  the  number  of  iterations  is  bounded  independent  of  h ,  H 

and  H/h.  (3.17b) 

•  The  number  of  iterations  for  the  multiplicative  Schwarz  method  is 
roughly  half  of  that  needed  for  the  additive  Schwarz  method.  (3.17c) 

The  local  linear  convergence  of  the  Algorithm  MSDM  is  shown  below  using  the  results 
in  [5]. 

Theorem  1  Let  pi™  are  the  n-th  step  MSDM  approximation  of  fia.  Assume  that 
Who.  —  hl\\vN  —  Hq  =  0,  •  •  • ,  n)  for  sufficiently  small  5.  Then,  we  have 

'  \K  -  A,lk.  <  ci(cy”|F(rf)  -  FW I.  (3.18) 

where  C\  is  a  positive  constant  and  C2  is  a  positive  constant  less  than  1.  (See  appendix 
for  the  proof  of  Theorem  1). 

3.4.  Computational  cost  of  the  Algorithms 

In  this  subsection,  the  computational  cost  of  the  proposed  Algorithms  will  be  estimated 
when  solving  (2.5)  and  (2.9)  for  the  forward  and  inverse  problem,  respectively.  We  will 
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call,  the  method  solving  (2.5)  and  (2.9)  without  any  domain  decomposition,  non-DDM 
compared  to  the  proposed  method.  The  system  matrices  for  each  equation  are  Nn  x  Nn 
and  Ne  x  Ne,  respectively.  Assume  that  we  require  O(Nf)  and  O(N'f)  floating  point 
operations  for  each  equation  with  non-DDM.  For  full  nonzero  matrix  we  must  take 
p  =  3  and  for  usual  non  sparse  matrix  we  can  take  2  <  p  <  3.  Suppose  that  we  use  d 
subdomains  with  equal  nodes  and  elements  for  the  proposed  algorithms. 

The  coarse  level  computation  needs  4~p  times  (in  two-dimension)  or  8~p  times  (in 
three-dimension)  less  computation  than  the  fine  level  computation  for  the  forward  and 
inverse  solvers.  Thus  if  we  neglect  the  coarse  level  computation,  the  computational 
cost  of  the  Algorithm  MODDM  is  MFdl~v  (for  one  computer)  and  MFd~p  (for  d 
parallel  computers)  times  that  of  non-DDM  method,  for  one  computer  and  for  d 
parallel  computers,  respectively.  It  is  well  known  that  the  number  of  iterations  MF 
is  bounded  independent  of  Nn,d  and  Nn/d  by  (3.17)  and  TMODDM  gives  a  good 
convergence  behavior  with  smaller  MF  than  MODDM.  Thus  we  have  chosen  MF  =  3 
when  TMODDM  is  used  in  the  inverse  solver  in  Section  4.2.  Thus,  by  using  multiple 
subdomains,  the  Algorithm  TMODDM  can  achieve  significant  decrease  in  computational 
requirements. 

The  total  computation  costs  for  non-DDM,  the  Algorithm  TMSDM  with  forward 
solver  the  Algorithm  TMODDM  in  case  of  one  computer,  d  parallel  computers,  and 
cP  parallel  computers  are  tabulated  in  Table  1.  When  we  use  d  parallel  computers, 
the  parallel  computing  is  applied  to  the  inverse  solver  only,  whereas  when  we  use 
d2  computers,  the  parallel  computing  is  applied  to  the  inverse  and  forward  solvers 
simultaneously.  Note  that  this  comparison  for  the  inverse  solver  is  based  on  the 
assumption  that  the  speed  of  data  communications  between  parallel  computers  is 
sufficiently  fast.  The  parallel  computation  is  not  treated  in  this  paper,  but  the  proposed 
algorithm  combined  with  the  parallel  computing  will  lead  to  much  more  efficient  results. 

4.  Numerical  Simulations 

In  this  section,  we  will  test  the  efficiency  of  the  two  algorithms:  Algorithm  TMODDM 
for  the  forward  problem  and  Algorithm  TMSDM  for  the  inverse  problem  using  simulated 
data. 

The  following  are  the  parameter  values  chosen  for  the  numerical  simulations. 

p!s  —  8 cm”1,  background  pa  =  0.05cm_1, 

tumor  pa  -  0.2 cm-1,  u  =  2n  *  100 MHz, 

c  =  3  *  10locm/sec,  a  =  1. 

In  Figure  1(a),  the  thin  lines  represent  20  x  20  fine  level  and  thick  lines  represent 
10  x  10  coarse  level  for  D  —  [0,6]  x  [0,6].  The  2x2  domain  decomposition  and  4x4 
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domain  decomposition  with  w  =  1  are  shown  in  Figure  1  (b)  and  (c).  3x1  domain 
decomposition  for  Q,  =  [0, 9, 6]  x  [0, 4.8]  with  24  x  18  node  points  is  shown  in  Figure  1 
(d). 

4-1.  Algorithm  TMODDM 

In  this  subsection,  we  consider  Algorithm  TMODDM  in  Figure  1,  Table  2,  and  Table 
3.  Let  0,  =  [0, 6]  x  [0, 6]  be  divided  into  20  x  20  square  elements.  The  distribution  of  pa 
is  shown  in  Figure  2(a),  where  the  white  region  is  the  background  with  pa  —  0.05cm"1 
and  the  black  region  is  the  anomaly  with  absorption  coefficient  pa  —  0.2 cm"1.  Let 
the  location  of  the  source  be  (6,3).  The  solution  of  (2.5)  by  LINPACK  is  in  Figure 
2(b).  The  L2(Q.)  errors  for  the  photon  density  between  LINPACK  solution  and  the 
solution  of  the  proposed  method,  are  displayed  in  Figures  2(c)-(e),  Table  2,  and  3.  In 
Table  2,  the  error  for  one-level  and  two-level  MODDM  is  compared  as  the  number  of 
iterations  increases  for  2  x  2  subdomains  as  in  Figure  1(b).  Although  TMODDM  4"p 
more  computation  than  MODDM,  the  convergence  behavior  of  the  proposed  method  for 
each  iteration  is  remarkably  good  as  compared  to  the  one-level  method.  In  Table  3,  the 
effect  of  the  number  of  subdomains  and  the  width  of  overlapping  region  is  tabulated. 
As  the  width  gets  larger,  the  error  decreases  with  a  few  exceptions.  This  verifies  the 
convergence  result  of  (3.17a).  In  Figure  2(c), (d),  and  (e),  the  L2  error  is  shown  after 
5  iterations,  after  the  coarse  level  correction  after  5  iterations,  and  after  6  iterations, 
respectively.  The  errors  for  each  step  are  2.51  x  10~9,  8.95  x  10"10,  and  8.31  x  10~ 11 , 
respectively. 

4-2.  Algorithm  TMSDM 

In  this  subsection,  we  evaluate  the  performance  of  the  algorithm  TMSDM  using 
simulated  data  and  Algorithm  TMODDM  as  a  forward  solver.  The  results  are  shown  in 
Figures  3-5.  The  reconstruction  of  diffusion  coefficient  k  for  known  pa  is  investigated  in 
[13].  In  this  paper,  the  reconstruction  of  pa  is  implemented  when  p's  is  given  as  in  [14]. 

Algorithm  TMSDM  is  implemented  for  various  parameters  including  type  of  domain 
Q,  pa,  source  and  detector  locations,  and  the  Tikhonov  regularization  parameter.  As 
a  forward  solver,  we  have  used  the  Algorithm  TMODDM  with  maximum  number  of 
iterations  set  to  3.  And  if  we  know  the  approximate  location  of  tumor  a  priori,  we 
only  updated  pa  only  on  some  small  region  of  interest  instead  of  all  subdomains  of  Q  in 
inverse  solver. 

The  maximum  number  of  iterations  for  the  coarse  level  correction,  the  subdomain 
correction,  and  non-DDM  is  10,  5,  and  10,  respectively.  As  shown  in  Table  1,  the 
computational  cost  of  the  proposed  method  is  less  than  that  of  non-DDM.  Figures 
??(a),  4(a),  and  5(a),  show  the  domain  fl  and  pa  images  used  in  simulations.  Figures 
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3(b),  4(b),  and  5(b),  the  location  of  the  sources  and  detectors,  ’o’  and  ’*’  mark  the 
location  of  sources  and  detectors,  respectively.  Figures  3(c)-(e),  4(c)-(f),  and  5(c)-(e), 
show  the  reconstructed  images  of  pa  for  various  parameters  and  the  L2{Cl)  norm  error 
between  the  original  pa  and  reconstructed  pa. 

In  Figure  3,  we  considered  Cl  =  [0, 6]  x  [0, 6]  domain  divided  into  400  pixels  with 
Nx  =  Ny  =  20.  We  choose  2x2  subdomains  of  the  same  size  with  overlapping  width 
1  as  shown  in  Figure  1(b).  Figure  3(a)  shows  the  original  image  of  pa.  20  sources 
are  interweaved  with  20  detectors  located  on  the  boundary  as  shown  in  Figure  3(b). 
The  image  is  reconstructed  using  the  Algorithm  TMSDM.  In  Figures  3(c)  and  3(d), 
the  reconstructed  images  are  shown  for  Tikhonov  regularization  parameter  a  =  0  and 
10~2,  respectively.  As  compared  to  the  image  with  no  regularization,  the  image  with 
a  =  10-2  appears  sharper  and  the  L2(Cl)  norm  error  is  decreased  from  0.145  to  0.129. 
The  reconstructed  image  for  the  non-DDM  after  10  iterations  is  shown  in  Figure  3(e). 
The  boundary  of  the  anomaly  does  not  appear  clearly  and  the  image  of  anomaly  is  larger 
than  the  original  image,  and  the  L2(Cl)  error  is  larger  than  that  of  the  DDM  approach 
shown  in  Figure  3(d). 

In  Figure  4,  the  same  domain  Cl,  pixels,  subdomains,  and  source  and  detector 
locations  are  considered  as  in  Figure  3.  But  the  anomaly  of  pa  is  located  at  the  upper 
right  quadrant  of  Cl  as  shown  in  Figure  4(a).  The  reconstructed  images  using  Algorithm 
TMSDM  with  a  =  0  and  a  —  10-2  are  shown  in  Figure  4(c)  and  Figure  4(d).  The 
reconstructed  image  when  pa  is  updated  only  in  4-th  subdomain  is  shown  in  Figure 
4(e).  The  reconstructed  image  for  the  non-DDM  is  shown  in  Figure  4(f).  The  L2(Cl) 
error  becomes  smaller  when  we  use  a  =  10-2  in  Figures  4  (d)-(f)  compared  to  a  =  0  in 
Figure  4  (c).  As  compared  to  Figure  4  (d),  when  we  only  update  the  4-th  subdomain,  the 
error  in  Figure  4  (e)  decreases  significantly.  Thus  if  a  priori  or  a  posteriori  information 
about  the  location  of  tumor  is  available,  it  is  possible  to  obtain  more  accurate  and  faster 
reconstruction  as  shown  in  Figure  4(e)  rather  than  Figure  4(d).  The  reconstructed  image 
for  the  non-DDM  in  Figure  4(f)  is  somewhat  blurred  and  the  boundary  is  not  so  clear 
as  compared  to  the  images  obtained  using  Algorithm  TMSDM. 

In  Figure  5,  we  considered  a  thin  slab  geometry  with  Cl  =  [0, 9.6]  x  [0, 4.8]  divided 
into  288  pixels  with  Nx  —  24  and  Ny  =  12.  We  have  choose  3x1  subdomains  of  similar 
size  with  overlapping  width  as  shown  in  Figure  1(d).  13  sources  on  the  bottom  and 
13  detectors  on  the  roof  of  Cl  are  located  as  shown  in  Figure  5(b).  The  reconstructed 
images  with  a  =  0  and  a  =  10-2  are  shown  in  Figure  5(c)  and  (d).  Figure  5(e)  is 
the  reconstructed  image  with  non-DDM.  The  reconstructed  image  with  regularization 
shown  in  Figure  5(d)  has  smaller  L2(Cl )  error  than  the  one  without  any  regularization 
shown  in  Figure  5(e).  As  in  Figure  3  and  4,  the  reconstructed  image  for  the  non-DDM 
shows  blurred  boundary  of  the  anomaly. 
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5.  Conclusion  and  discussion 

We  presented  Algorithm  TMODDM  and  Algorithm  TMSDM  for  the  forward  and  the 
inverse  problem  solution  of  DOT.  The  algorithms  of  these  two  methods  are  explained 
and  their  convergences  are  shown.  Numerical  implementations  of  the  algorithms  are 
presented  for  various  parameters  such  as  geometry,  source-detector  location,  Tikhonov 
parameter,  the  number  and  width  of  subdomain.  For  Algorithm  TMODDM,  the  forward 
solver,  as  overlapping  width  grows  the  convergence  becomes  better  for  the  4  and  16 
subdomains,  and  the  iteration  number  for  the  proposed  method  decreases  fundamentally 
compared  to  the  one-level  method.  For  Algorithm  TMSDM  combined  with  Algorithm 
TMODDM,  Tikhonov  parameter  1CT2  produces  less  L2(D)  error  for  the  absorption 
coefficient  than  other  Tikhonov  parameters.  The  images  of  the  coefficient  by  non- 
DDM  are  blurred  and  makes  higher  L2(Q)  error  than  the  proposed  method.  Using 
Algorithm  TMSDM  in  a  region  of  interest  produces  better  quality  images  and  less  error 
as  compared  to  the  Algorithm  TMSDM  applied  in  all  subdomains. 

The  computational  costs  for  the  two-level  methods  are  compared  among  non-DDM, 
the  proposed  method  in  1,  d,  and  d?  parallel  computers,  where  d  is  the  number  of 
subdomains.  Thus,  with  the  aid  of  parallel  computing  for  d2  computers,  the  proposed 
approach  speeds  up  the  computation  with  d~p  times  as  compared  to  the  non-DDM, 
where  usually  p  is  between  2  and  3  depending  on  the  system  matrix  size  and  matrix 
solver.  If  a  region  of  interest  is  determined  in  some  subdomain  by  a  priori  information 
obtained  from  secondary  imaging  modality  such  as  X-ray,  CT,  MRI  or  a  posteriori 
information  from  a  coarse  level  image,  Algorithm  TMSDM  can  be  implemented  in  the 
region  of  interest  resulting  in  a  d~v  speed-up  of  computation  only  in  d  computers. 

6.  Appendix:  the  proof  of  Theorem  1 

More  detailed  proof  of  the  Theorem  in  general  case  including  the  concrete  definitions  of 
Jacobian  and  Hessian,  and  their  error  estimates,  will  be  studied  in  [15].  In  this  appendix, 
we  use  absolute  intensity  r  =  as  a  measurement  data  for  real  valued  $  instead  of 
log  intensity  in  (2.2),  and  assume  that  k  =  ^r,  which  means  =  0. 

For  i-th  source  ri)  j- th  detector  r,-  and  elements  Tk  and  TJ,  let  btj  =  Tiffpa  +  8pa)  — 
Tij(jia).  The  Jacobian  matrix  J  and  Hessian  tensor  IT  at  a  given  pa  are  given  by 

J(fxa)[5pa\(ij,k)  =  [  G(ri,r,)G(r',rJ)J/xa(r,)d^ 

Jrk 

H{pa){Spa,dHa](ij,k,l) 
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where  G(-,rf)  is  the  Green  function  for  pa  and  a  point  source  rj.  Define 

||J||  :=m&x\J(ij,k)\  <  sup||G(-,r)||*a(n)  ||^a||LOO(n) , 

ren 

Ill'll  :=  max \H(ij,k,l)\  <  sup  ||G(-,r)||®a(n)  ||<W£oo(n) . 

reel 

Lemma  2  Assume  that  5fj,a  is  small  enough  to  satisfying  \\H\\  ||  J\\  ||^a||LO0^^  <  | 
and  suprgQ  G(-,r)  <  §supren  ||G(-,r)||£2(n),  where  G(-,r)  is  the  solution  of  (2.1) 

for  the  absorption  coefficient  pa  +  8 pa.  Then  F"  satisfies 

§  IIWIv*.  <  (8^yF"(x)(8^a).  (6.19) 

Proof  F'  and  F"  are  computed  as  follows: 

F'iPc)  =  +  aji,,. 

F"M  =  Jfjt.YJtp.)  +  HMbM  +  al. 

By  Born  expansion 

K\  =  \G(ri+Ns,rj)  -  G{ri+Na,rj) \  =  |  /  G(rhr')G(r' ,rj)8ha\, 

Jn 

where  G  is  the  solution  of  (2.1)  for  fsa  +  8pa.  Using  the  assumption  and  the  definition 
of  J,  we  get  \bij\  <  § 1|  J\\  pMa||ioo(n).  Therefore  we  get 

§IIWIy„.  <  (o- |  ||fl(M.)||  II JWII  11^11^(0,)  ll^ll2^ 

<  (6.20) 

Thus  we  have  proved  the  lemma. 

The  inequality  (6.19)  is  called  local  strong  convexity.  With  this  property  and 
the  theorem  in  [5],  we  prove  Theorem  1  for  6  is  sufficiently  small  such  that  5pa  = 
iff  —  mua,  q  =  0,  •  •  • ,  n  satisfies  the  assumption  of  Lemma  2. 

Proof  of  Theorem  1  The  proof  of  Theorem  1  is  the  modification  of  the  proof  in  [5]. 
By  (3.14), 

vNe  =  vl  +  ---  +  v«e  =  w1Ne  +  ---  +  w«e 

and  Wffe,p  =  are  mutually  disjoint  and  Wfe  C  Vfe.  And  let  and  ii,an+k^d 

be  the  minimizing  solution  of  TMSDM  at  n-th  step  and  p-th  subdomain  0“ .  Define 
zf  =  (ha  ~  han+*)\np  e  WpNe  and  en+p/d  =  pan+d  -  pian+2ir  e  V£e.  Then  we  get 

d 


P= 1 


ha  ha  II  ha  ha  ||yNe  — 
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and 

en+d  =  arg  min Vp€yNepF(pan+E^  +  Vp). 

(6.22)  implies 

(F'(pan+E?  +  en+d),vp  -  en+ty  >  0  for  all  vp  G  VNf. 

Using  (6.19)  and  Taylor  expansion  for  F  we  get 


a 


F(w)  —  F(v)  >  { F'(v ),  w  —  v)  +  —  ||iu  —  u|| yN 
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(6.22) 

(6.23) 

(6.24) 


for  all  v,w  G  Vmb  such  that  ||u  —  w\\Vn  <  S.  Taking  vp  =  0  in  (6.23)  and  inserting 
w  =  pan+ ^  and  v  =  pan+*  we  get 


F(pan+^)  -  F(pan+Zd)  > - 


en+z 


>0 


vNe 


and 


FW)  -  f(M«"+1)  >  E  iF(^n+C?)  -  F(, M/+5)) 

P=  1 

^4i< 


|e"+§ 


>  0. 


VnF 


Inserting  w  =  pan  and  v  =  /ia  into  (6.24),  we  get 


a 


F(fian)  -  FW  =  W,  tf."  ~  Ha)  +  J  |k"  ~ 


a 


>  J  Ik"  -  >  0. 

Finally  we  get 

o  <  f(ft,”+1)  -  F(/i.)  <  <F'(^"+1),^”+1  -  (using  (6.24)) 
=  (  F'(X.”+1),  E  e”+5  +  /*.”  -  A.)  =  E  (F'(fc"+1,  s"+5  -  jj 


p=l 


p=l 


a 

<  £  (F'(pan+1)  ~  F'(/xan+*),en+*  -  zp")  (using  (6.23)) 

p=i 

=  E  E  (irV."+5)  -  F'(ft.”+T),e“+  5  - 

P=1  i=p+l 

=  E  E  (F"(»r)e"+ V+!  -  *”) 


p=i  i=p+i 

for  some  y"  which  lies  between  pan+^  and  pan+L^ 


(6.25) 


(6.26) 
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d  „  \  1/2  /  d  .  \  V2 


< 


a 


E 

j= 1 


en+~d 


E 

vP=l 


on+ 


d  —  z 


(using  (6.19)) 


1/2' 


<  2 (F(pan)  -  F(pan+'))  +  2 y/FM  ~  F(pan+1)\/F(pan)  -  F(pa) 
(using  (6.21),  (6.25),  and  (6.26)  ). 


Let  dn  =  F(pan)  —  F{pa),  then  the  above  equation  and  Cauchy-Schwarz  inequality 
implies 


d 


n+1 


E  2 {dn  dn-)_i)  T  2  dn  dn-!  i  \j dn 

E  (2  +  t  )(dn  —  dn+i)  +  2pdn, 


for  all  n  >  0.  Thus  we  get 

dn+ i  <  4/i1 2 3 4  +  4/i  +  1 
dn  ~  6//  +  1 


(6.27) 


The  left  hand  side  of  (6.27)  takes  minimum  value  §  <  1  when  p  =  \-  Thus  if  we  take 
C\  =  ~  and  C'i  =  Theorem  1  is  proved. 
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Method 

Computational  cost 

non-DDM 

0(A®  +  0(JV?) 

TMSDM  on  1  computer 

MFdl~pO{Np)  +  d1~pO(Np) 

TMSDM  on  d  computers 

MFdl'pO{Np)  +  d~pO(Np) 

TMSDM  on  d2  computers 

MFd~pO(Np)  +  d~pO(Np) 

Table  1.  The  comparison  of  the  computational  costs  for  the  non-DDM  method, 
TMSDM  for  one,  d,  and  d2  computers. 


Method\  Iteration 

1 

2 

3 

4 

5 

6 

one-level 

.318e-01 

.148e-02 

.555e-03  .143e-03 

.235e-04 

.254e-05 

two-level 

.236e-01 

.466e-03 

.717e-05 

.104e-06 

.251e-08 

.831e-10 

Table  2.  The  comparison  of  L2  error  with  respect  to  the  LINPACK  solution,  which 
is  shown  in  Figure  2(b),  for  one-level  and  two-level  MODDM,  as  iteration  increases. 


Subdomains  \  Iteration 

1 

2 

3 

4 

5 

6 

2x2  subdomains,  width  1 

.236e-01 

.466e-03 

.717e-05 

.104e-06 

.251e-08 

.831e-10 

2x2  subdomains,  width  2 

.248e-02 

.306e-05 

.553e-07 

.204e-09 

.478e-12 

.529e-14 

2x2  subdomains,  width  3 

.195e-03 

.218e-06 

.101e-09 

.571e-13 

,167e-15 

.163e-15 

4x4  subdomains,  width  1 

.237e-01 

.482e-03 

.808e-05 

.208e-06 

.894e-08 

.185e-09 

4x4  subdomains,  width  2 

.329e-02 

.106e-04 

.503e-06 

.209e-08 

.351e-10 

,327e-12 

4x4  subdomains,  width  3 

.706e-03 

.231e-05 

.440e-08 

.136e-10 

,107e-13 

.181e-15 

Table  3.  The  effect  of  multiple  subdomains  and  overlapping  width  to  the  convergence. 
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(a) 


(b) 


(c)  (d) 

Figure  1.  (a)  The  lattice  of  thin  lines  is  20  x  20  mesh  of  ft  =  [0, 6]  x  [0, 6],  and  thick 
lines  represents  the  coarse  level  of  0,  (b)  2x2  domain  decomposition  of  ft  =  [0, 6]  x  [0, 6] 
with  20  x  20  mesh,  (c)  4  x  4  domain  decomposition  of  ft  =  [0, 6]  x  [0, 6]  with  20  x  20 
mesh,  (d)  3  x  1  domain  decomposition  of  ft  =  [0, 9.6]  x  [0, 4.8]  with  24  x  12  mesh. 
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Figure  2.  (a)  is  the  image  of  /j.a.  The  white  region  represent  background  tissue 
with  value  0.05cm”1  and  the  black  region  represent  anomaly  which  have  0.2cm”1 
absorption  coefficient.  The  photon  density  with  point  source  at  (6,3)  solved  by 
LINPACK  is  depicted  in  (b).  The  error  between  LINPACK  solution  and  the  result 
after  5th  iteration,  coarse  level  after  5th  iteration,  and  6th  iteration  attained  from 
Algorithm  TMODDM  is  presented  in  (c),  (d),  and  (e),  respectively. 


Figure  3.  The  reconstruction  of  absorption  coefficient  using  Algorithm  TMSDM  on 
Cl  =  [0,6]  x  [0,6]  with  2x2  subdomains,  20  detectors,  20  sources,  and  400  pixels 
(a)The  image  of  fj.a  to  be  reconstructed.  White  region  represent  background  tissue 
with  fia  =  0.05cm-1  and  black  region  represent  anomaly  with  jia  =  0.2cm-1  (b)The 
location  of  20  sources  and  20  detectors  (c)The  reconstructed  image  of  the  absorption 
coefficient  with  no  regularization  using  Algorithm  TMSDM  resulting  in  L 2  norm  error 
0.145  (d)The  reconstructed  image  of  the  absorption  coefficient  with  a  =  10-2  using 
Algorithm  TMSDM  resulting  in  L2  norm  error  0.129  (e)The  reconstructed  image  of 
the  absorption  coefficient  with  a  =  10-2  for  non-DDM  resulting  in  L2  norm  error  0.142 
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Figure  4.  The  reconstruction  of  the  absorption  coefficient  on  fi  =  [0, 6]  x  [0, 6]  with 
2x2  subdomains,  20  detectors,  20  sources,  and  400  pixels,  (a)  The  image  of  ga  to 
be  reconstructed.  White  region  represent  background  tissue  with  ya  =  0.05cm_1  and 
black  region  represent  anomaly  with  pa  —  0.2cm-1.  (b)  The  location  of  20  sources  and 
20  detectors  (c)  The  reconstructed  image  of  the  absorption  coefficient  using  Algorithm 
TMSDM  with  a  =  0  resulting  in  L 2  norm  error  0.0949  (d)  The  reconstructed  image 
of  the  absorption  coefficient  using  Algorithm  TMSDM  with  a  =  10-2  resulting  in 
L2  norm  error  0.0824  (e)  The  reconstructed  image  of  the  absorption  coefficient  using 
Algorithm  TMSDM  only  on  fourth  subdomain  with  a  —  10-2  resulting  in  L2  norm 
error  0.0755  (f)  The  reconstructed  image  of  the  absorption  coefficient  for  non-DDM 
with  a  =  10-2  resulting  in  L 2  norm  error  0.0926 


24 


Two-level  Domain  Decomposition  Methods  for  Diffuse  Optical  Tomography 


(e) 


Figure  5.  The  reconstruction  of  absorption  coefficient  on  fl  =  [0, 9.6]  x  [0, 4.8]  with 
3x1  subdomains,  13  detectors,  13  sources,  and  288  pixels,  (a)  The  image  of  Ha  to 
be  reconstructed.  White  region  represent  background  tissue  with  ya  =  0.05cro-1  and 
black  region  represent  anomaly  with  \ia  =  0.2cm_1.  (b)  The  location  of  13  sources  and 
13  detectors  (c)  The  reconstructed  image  of  the  absorption  coefficient  using  Algorithm 
TMSDM  with  a  =  0  resulting  in  L2  norm  error  0.156.  (d)  The  reconstructed  image 
of  the  absorption  coefficient  using  Algorithm  TMSDM  with  a  —  10-2  resulting  in 
L2  norm  error  0.138.  (e)  The  reconstructed  image  of  the  absorption  coefficient  for 
non-DDM  with  a  —  10-2  resulting  in  L2  norm  error  0.141. 
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ABSTRACT 

In  this  paper,  we  explore  domain  decomposition  algorithms  for 
the  inverse  DOT  problem  in  order  to  reduce  the  computational 
complexity  and  accelerate  the  convergence  of  the  optical  image 
reconstruction.  We  propose  a  combination  of  a  two-level  multi¬ 
grid  algorithm  with  a  modified  multiplicative  Schwarz  algorithm, 
where  a  conjugate  gradient  is  used  as  an  accelerator  to  solve  each 
sub-problem  formulated  on  each  of  the  partitioned  sub-domains. 
For  our  experiments,  simulated  phantom  configuration  with  two 
rectangular  inclusions  is  used  as  a  testbed  to  measure  the  compu¬ 
tational  efficiency  of  our  algorithms.  No  a  priori  information  about 
the  configuration  is  assumed  except  for  the  source  and  detector  lo¬ 
cations.  For  the  application  of  our  modified  Schwarz  algorithm 
alone,  we  observe  an  increase  in  efficiency  of  100%  as  compared 
to  the  conjugate  gradient  solution  obtained  for  the  full  domain. 
With  the  addition  of  the  coarse  grid,  this  efficiency  rises  to  400%. 
The  coarse  grid  also  serves  to  improve  the  overall  appearance  of 
the  reconstructed  image  at  the  boundaries  of  the  inclusions. 

1.  INTRODUCTION 

Inverse  diffuse  optical  tomography  (DOT)  problem  involves  es¬ 
timation  of  the  optical  properties  of  biological  tissues,  which  are 
pertinent  to  tissue’s  physiological  and  biochemical  state.  The  most 
prominent  applications  of  DOT  are  in  detecting  tumors  in  the  breast, 
monitoring  brain  activity,  and  detecting  brain  tumors  and  hemor¬ 
rhages. 

DOT  poses  a  computationally  challenging  inverse  problem. 
Thereby,  realization  of  real-time  diffuse  optical  tomographic  imag¬ 
ing  requires  computationally  viable  reconstruction  algorithms  that 
provide  accurate  quantitative  results.  In  this  work,  we  address  the 
computational  complexity  of  the  inverse  problem  by  proposing  a 
two-level  domain  decomposition  procedure.  Domain  decomposi¬ 
tion  methods  convert  the  inverse  problem  into  smaller-size  prob¬ 
lems  that  are  easier  to  handle.  Herein  we  apply  a  modified  alter¬ 
nating  Schwarz  method  with  conjugate  gradient  (CG)  algorithm 
as  the  iterative  solver  to  accelerate  the  solution.  We  then  extend 
the  uni-level  problem  formulation  to  a  2-level  problem  to  include 
a  coarse  grid  correction  in  an  attempt  to  improve  quantitative  ac¬ 
curacy  and  further  accelerate  the  convergence. 

The  Schwarz  algorithm  begins  by  partitioning  the  domain  into 
two  or  more  overlapping  sub-domains.  The  problem  is  then  di¬ 
vided  into  subproblems  on  each  of  these  sub-domains.  Multiplica¬ 
tive  Schwarz  algorithm  successively  solves  the  localized  problems 
in  each  sub-domain,  every  time  using  the  latest  solution  available 
to  initialize  the  next  sub-problem.  Therefore,  at  every  iteration 


one  sub-problem,  initialized  with  the  solution  from  the  last  itera¬ 
tion,  is  solved,  whose  solution  is  then  used  to  initialize  the  next 
sub-domain  and  so  on.  We  follow  the  multiplicative  Schwarz  al¬ 
gorithm  using  a  fixed  number  of  CG  iteration  to  approximate  the 
solution  for  each  sub-domain.  For  full  explanation  of  Schwarz  al¬ 
gorithms  see  [1, 2]. 

Single  level  Schwarz  algorithms  are  not  well  suited  for  prob¬ 
lems  exhibiting  low  frequency  errors  and  suffer  from  low  conver¬ 
gence  rate.  Multi-grid  methods,  on  the  other  hand,  use  a  hierarchy 
of  grids  at  different  scales  to  accelerate  the  convergence  of  stan¬ 
dard  iterative  methods.  The  fundamental  idea  behind  all  multi-grid 
methods  is  to  combine  computations  done  on  different  grid  scales 
in  order  to  eliminate  the  error  components  of  the  finest  grid,  where 
the  original  problem  has  been  formulated.  This  is  achieved  by  ap¬ 
proximating  the  smoothed  fine  grid  error  on  a  coarser  grid  where 
it  can  be  accurately  represented.  The  approximated  error  on  the 
coarse  grid  appears  to  be  more  oscillatory,  hence  can  be  further 
eliminated  by  the  iterations  on  the  coarse  grid.  The  solution  ob¬ 
tained  for  the  error  on  the  coarse  grid  is  then  interpolated  to  the 
fine  grid  to  correct  the  current  fine  grid  solution  estimate.  There 
is  no  unique  way  of  formulating  an  inverse  problem  in  a  multigrid 
framework.  The  formulation  needs  to  be  specific  to  the  problem  at 
hand,  that  is  it  must  be  tuned  to  address  the  requirements  asserted 
in  the  problem.  Note  that,  carefully  designed  multi-grid  solvers 
have  the  potential  of  solving  inverse  problems  with  N  unknowns 
within  O(N)  work  load  [3],  which  makes  them  the  most  efficient 
solvers  for  many  kinds  of  mathematical  problems. 

Multigrid  methods  have  been  applied  for  the  DOT  problem 
in  the  past,  where  the  inverse  problem  is  formulated  on  a  hier¬ 
archy  of  regular  rectangular  grids  [4].  Domain  decomposition 
has  also  been  previously  proposed  for  Bayesian  formulation  of  in¬ 
verse  DOT  problem  [5],  However,  the  nature  of  the  decomposition 
was  different  in  that  work,  where  non-overlapping,  and  hence  non- 
Schwarz  type,  ’’sliding  window”  decomposition  was  used.  In  our 
previous  work,  we  have  proposed  Fast  Adaptive  Composite-grid 
(FAC)  algorithms  for  Region-of-Interest  DOT  [6, 7],  which  implic¬ 
itly  pursued  a  domain  decomposition  with  the  aid  of  a  priori  infor¬ 
mation  about  the  medium  of  interest.  Note  that  FAC  can  be  viewed 
as  a  Schwarz-like  domain  decomposition  method  in  terms  of  fully 
overlapping  sub-domains,  hence  achieving  fast  convergence  with 
low  cost  by  the  use  of  coarse  grid  with  substantially  fewer  points 
in  the  overlap  region  [8], 

In  this  work,  we  assume  that  no  a  priori  information  is  avail¬ 
able  about  the  unknown  image.  We  have  used  the  location  of 
sources  and  detectors  to  determine  the  overlapping  sub-regions. 
We  formulate  the  inverse  problem  on  the  fine  grid  as  two  smaller 


size  inverse  problems  on  overlapping  sub-regions.  CG  accelerator 
is  used  to  approximate  the  solution  on  one  sub-region.  This  solu¬ 
tion  is  used  to  initialize  the  optimization  in  the  other  sub-region. 
The  fine  grid  iterations  are  followed  by  a  coarse  grid  correction 
scheme,  where  the  inverse  problem  aims  to  solve  the  residual  equa¬ 
tion  formulated  on  the  global  coarse  grid,  which  is  of  relatively 
smaller  size.  A  number  of  2-grid  cycles  are  run  to  further  improve 
the  accuracy  of  the  reconstruction. 

2.  PROBLEM  FORMULATION 

2.1.  Forward  Model 

The  forward  model  for  DOT  is  based  on  simplifying  assumptions 
applied  to  radiative  transport  equation  which  results  in  a  form  of 
photon  diffusion  equation.  In  frequency  domain,  the  diffusion 
equation  is  given  by: 

-DV2$(r)  -  3—$(r)  +  pa{r)$(r)  =  S  (1) 

c 

where  D  is  the  diffusion  coefficient,  c  is  the  speed  of  light  and 
Pair)  on  the  entire  domain  Cl  C  R2  is  the  spatially  varying  ab¬ 
sorption  coefficient.  S  stands  for  the  point  source  located  at  r  = 
rs.  In  this  work,  we  focus  on  the  reconstruction  of  absorption 
coefficients  pa  of  the  medium,  hence  we  assume  that  the  diffu¬ 
sion  coefficient  D  has  a  spatially  uniform  distribution.  We  have 
employed  the  perturbation  approach  [9]  around  a  spatially  invari¬ 
ant  optical  background  with  a  first  order  Rytov  approximation  to 
solve  the  inverse  problem.  The  cell-centered  discretization  on  the 
grid  Clh,  yields  a  system  of  linear  equations  relating  the  differential 
absorption  coefficients  Spain)  to  the  measurements: 
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where  y(j  denote  the  real  part  of  the  measurement  at  ith  source 
and  jth  detector  at  frequency  /*,,  W/jkh  is  the  weight  for  the  Ith 
pixel  for  ij  source-detector  pair,  and  Sfia(ri)  is  the  differential  ab¬ 
sorption  coefficient  for  lth  pixel.  We  can  denote  this  model  system 
succinctly  as: 

y  =  Whxh  (3) 

where  Wh  is  the  weight  matrix  and  y  €  RM,  x  e  R'v,  Wh  € 
VM  xN .  N  is  the  number  of  grid  points  on  Clh  and  M  is  the  total 
number  of  measurements. 

2.2.  Inverse  Problem 

We  formulate  the  discrete  inverse  problem  to  yield  a  minimum 
least  squares  solution  for  the  differential  absorption  coefficients 
xh  on  Clh. 

Xls  =  argmin  4'(xh)  =  argmin  ||y  —  Whxh ||2  (4) 

Xh  Xh 

where  1 1  •  1 1  denotes  the  Euclidean  norm. 


2.2. 1.  Uni-Level  Domain  Decomposition  Algorithm 


Fig.  1.  Decomposition  of  Cl 

In  this  work,  we  decompose  the  domain  Cl  into  two  overlap¬ 
ping  sub-domains  fli  and  CI2  such  that  Cl  =  Hi  U  fi2  and  Go  is 
the  overlapping  region,  that  is  Do  =  Dj  n  Cl2  (Figure  1).  We  dis¬ 
card  those  source-detector  pairings  that  would  produce  a  coupling 
between  the  two  sub-domains  outside  the  overlapping  region.  In 
other  words  we  remove  any  measurement  that  is  due  to  a  source  in 
f2i\fi2  and  to  a  detector  in  O2  \Oi  or  vice  versa.  As  a  result,  the 
measurement  vector  y  e  RM  for  M  <  M  becomes 

\ 

VnQ  (5) 

3/n2\ni  / 

Regrouping  the  measurements  in  two  vectors  yi  and  3/2  yields 
yi  =  (j/nAnJ  m 0)T  and  2/2  =  (yn0  |  l/n2\n1)T-  Similarly  the 
differential  absorption  coefficients  on  the  sub-domains  are  grouped 
as  Xi  and  x2,  such  that  x2  is  a  finite  dimensional  approximation 
of  differential  absorption  coefficients  on  fif  and  x2  on  Cl2.  Asa 
result  we  can  formulate  two  sub-domain  problems  as  follows 

Xi,ls  =  argmin  ||t/i  -  WiXif 

*1 

&2,ls  =  arg  min  \\y2  -  W£x 2  ||2  (6) 

*2 

Instead  of  estimating  xh  using  the  formulation  given  in  equation  4, 
we  propose  a  procedure  that  follows  minimization  of  the  two  ob¬ 
jective  functionals  formulated  on  the  two  sub-domains,  with  re¬ 
duced  number  of  measurements.  This  results  in  significant  reduc¬ 
tion  in  the  size  of  the  inverse  problem  and  consequently  the  com¬ 
putational  complexity  of  the  overall  inverse  DOT  problem. 

A  conjugate  gradient  (CG)  algorithm  is  utilized  to  accelerate 
the  solution  on  each  sub-grid.  The  solution  update  obtained  after  a 
sweep  of  CG  iterations  on  one  domain  is  restricted  to  the  overlap¬ 
ping  region  on  the  grid  f2o  by  the  discrete  operator  I :  fif  —>  CIq. 
The  restricted  estimates  are  then  used  to  update  the  x2  estimates 
on  Gq  ,  which  is  followed  by  iterations  on  to  yield  a  solution  on 
entire  discrete  domain  Cl2.  A  v  number  of  cycles  is  applied  until  a 
desired  level  of  convergence  is  achieved. 

This  approach  takes  advantage  of  the  reduced  size  of  inverse 
problem  formulation  by  decomposing  the  domain  and  the  associ¬ 
ated  measurements.  Initialization  in  one  sub-problem  solution  by 
the  current  estimates  in  the  other  one  and  use  of  CG  algorithm  are 
other  important  factors  facilitating  the  acceleration  of  the  solution. 

2.2.2.  Integrating  Coarse  Grid  Correction 

Single  level  Schwarz  algorithms  are  not  well  suited  for  problems 
exhibiting  low  frequency  errors.  CG  iterations  tend  to  smooth 


the  error  in  the  solution.  Even  though  CG  algorithm  accelerates 
the  convergence  to  the  solution  xh  on  flh,  it  is  unable  to  further 
eliminate  the  error  with  smooth  content.  This  necessitates  the  use 
of  a  coarse  grid  correction  scheme,  which  enables  elimination  of 
smooth  error  by  restricting  it  onto  a  coarser  grid  and  correcting  the 
fine  grid  solution  by  interpolating  the  error  in  the  coarse  grid. 

After  the  end  of  CG  cycles  on  the  sub-problems,  the  current 
estimates  Xi  and  are  concatenated  to  yield  the  overall  solution 
estimate  xh  on  The  error  between  the  actual  solution  x*  and 
the  current  solution  estimate  x h  is  given  by  eh  =  x*  -  xh.  As¬ 
suming  that  the  error  eh  on  Qh  is  smoothed  well  enough  by  CG 
iterations,  we  can  write  a  coarse  grid  approximation  for  this  error, 
such  that 


h  rh  2h 

e  —  h  h.e 


where  l£h  is  the  interpolation  operator.  We  have  selected  I%h  : 

RN/4  RN  as  4(/2h)T  where  j2h  .  RiV  - >  RN/ 4  ^  (he  fuU 

weighting  operator  in  2D  [10].  As  a  result,  the  objective  functional 
given  in  equation  4  can  be  re-written  as 


\\y-Whxh\\2 


\\y-Wh{xh  +  eh)\\2 
\\y  -  Wh(xh  +  I$he2h)f 
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which  is  the  coarse  grid  objective  functional  defined  on  the  coarse 
grid  fl2h  with  grid  size  of  2h.  r  =  y  -  Whxh  is  called  the  resid¬ 
ual,  e2h  is  the  coarse  grid  error,  and  W2h  =  Whl2h  is  the  coarse 
grid  operator.  CG  iterations  on  the  coarse  grid  error  e2h,  which 
is  initially  assigned  to  zero,  eliminate  the  high-frequency  compo¬ 
nents  of  e2h  which  appeared  to  be  smooth  on  Qh.  This  results  from 
the  fact  that  low  frequency  components  appear  more  oscillatory  on 
coarse  grids  as  compared  to  fine  grids. 

The  error  estimate  e2h  on  fi2h  can  then  be  interpolated  to  the 
fine  grid  to  correct  the  fine  grid  solution  estimate  xh 

xh  4-  xh  +  l£he2h  (8) 


A  predefined  number  of  2-grid  cycles  are  run  to  further  improve 
the  solution  accuracy.  A  pseudo-code  of  the  overall  2-level  domain 
decomposition  algorithm  is  given  below. 


Algorithm  1  Two-Grid  Domain  Decomposition  Algorithm 

1:  Wi,W2,yi,y2*-partition(Wh,y) 

2:  W2h  <—  restrict{Wh)  {We  only  need  to  generate  U2h  once} 
3:  xh  <—  initialize(xh) 

4:  repeat 

5:  Xi  <—  CG(Wi,yi,xi)  {Smoothing  on  fli } 

6:  *2  <—  CG(W2,y2,X2)  {Smoothing  on  fl2} 

7:  r  <—  y  —  Wh x h  {Calculate  residual.} 

8:  e2h  <—  CG(W2h,  r,  e2h)  {Solve  for  coarse  grid  error} 

9:  eh  <—  interpolate(e2h ) 

10:  xh  <-  xh  +  eh 
11:  until  convergence 


3.  EXPERIMENTAL  RESULTS 

We  performed  test  cases  against  a  simple  simulated  phantom  con¬ 
figuration  consisting  of  two  rectangular  inclusions  as  seen  in  Fig. 


3.  The  resolution  of  our  test  image  was  20  pixels  by  40  pixels  and 
the  number  of  sources  and  detectors  used  was  17  and  33  respec¬ 
tively  resulting  in  a  561  by  800  weight  matrix  Wh.  We  compared 
four  test  cases:  (1)  least  squares  solution  using  conjugate  gradient 
on  flh  for  the  full  inverse  problem,  (2)  least  squares  solution  us¬ 
ing  conjugate  gradient  on  flh  with  the  reduced  source  and  detector 
configurations,  (3)  uni-level  domain  decomposition  on  the  reduced 
source  and  detector  configuration,  and  (4)  two-level  domain  de¬ 
composition  with  modified  multiplicative  smoother  on  the  reduced 
source  and  detector  configurations.  By  discarding  those  measure¬ 
ments  that  couple  source  on  one  sub-domain  and  the  detector  on 
the  other,  we  effectively  reduce  the  dimension  of  the  weight  matrix 
Wh. 

All  algorithms  were  implemented  in  MATLAB.  The  results 
were  compared  using  the  same  conjugate  gradient  and  decompo¬ 
sition  codes.  Figure  2  shows  a  plot  of  square  error  between  the 
actual  and  estimated  image  versus  floating  point  operations  (flops) 
required  for  reconstruction.  The  square  error  was  calculated  by 
taking  the  pixel  by  pixel  difference  between  the  true  image  (see 
Figure  3)  and  the  reconstructed  image  then  taking  the  sum  of  the 
squares  of  those  differences.  As  shown  in  Figure  2,  there  is  up  to 
an  average  of  100%  increase  in  efficiency  using  the  uni-level  do¬ 
main  decomposition  algorithm  compared  with  case  (1).  With  the 
addition  of  coarse  grid  correction,  there  is  approximately  400% 
increase  in  efficiency.  We  observed  that  the  error  curve  tends  to 
settle  faster  for  single  level  methods.  This  is  to  be  expected  since 
the  motivation  for  coarse  grid  correction  is  to  compensate  for  the 
smoother’s  inability  to  properly  handle  globalized  (i.e.  low  fre¬ 
quency)  errors. 


Fig.  2.  Square  Error  vs.  No.  of  Flops 
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Fig.  3.  Simulated  Phantom  Configuration 

Figures  4, 5, 6,  and  7  shows  the  reconstructed  images  for  each 
method  after  1200  iterations.  As  can  be  inferred  from  Figures  5 
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Fig.  4.  Direct  Method:  1200  Iterations 


Fig.  6.  One  Grid  D.D.:  1200  Iterations 


Fig.  5.  Direct  Method  w/  Reduced  S-D  Pairs:  1200  Iterations 


Fig.  7.  Two  Grid  D.D.:  1200  Iterations 


and  6,  the  uni-level  domain  decomposition  method  and  the  direct 
CG  on  the  reduced  source-detector  configurations  had  nearly  equal 
error  after  1200  iterations,  but  the  computational  time  was  cut  by 
half  for  the  domain  decomposition  case.  The  effect  of  coarse  grid 
correction  on  the  domain  decomposition  algorithm  is  clear  from 
Figure  7.  By  incorporating  coarse  grid  correction,  the  algorithm 
is  able  to  handle  the  smooth  errors  around  the  boundaries  of  the 
inclusions.  In  effect,  we  end  up  with  a  much  sharper  picture  and 
better  results  around  these  boundaries. 

4.  CONCLUSION 

In  this  paper,  we  investigated  the  effectiveness,  in  terms  of  com¬ 
putational  efficiency,  of  applying  Schwarz  type  domain  decom¬ 
position  method  in  the  solution  of  linearized  DOT  inverse  prob¬ 
lem.  We  modified  the  classic  Schwarz  algorithm  where  CG  was 
used  to  accelerate  the  convergence  to  the  least  squares  solution 
for  each  sub-problem  in  the  place  of  the  usual  Gauss-Siedel  itera¬ 
tions.  In  our  two  grid  domain  decomposition  algorithm,  CG  was 
applied  to  each  sub-domain  for  several  iterations  (to  get  a  good 
smoothing  effect).  The  residual  from  the  entire  domain  was  then 
restricted  to  the  coarse  grid  to  find  the  coarse  grid  error  which  was 
then  interpolated  back  to  the  fine  grid  for  error  correction.  The 
addition  of  the  domain  decomposition  reduced  the  computational 
burden  while  addition  of  the  coarse  grid  allowed  for  correction 
of  global  (low  frequency)  errors.  In  terms  of  computational  effi¬ 
ciency  as  compared  with  CG  solution  for  the  full  problem,  uni¬ 
level  domain  decomposition  method  showed  100%  increase  while 
two-level  method  saw  nearly  400%.  As  a  result  of  these  prelim¬ 
inary  findings,  we  feel  that  domain  decomposition  coupled  with 
multi-grid  methods  is  a  viable  option  for  decreasing  the  computa¬ 
tion  complexity  of  the  DOT  inverse  problem  even  when  there  is  no 
parallelization  of  the  problem.  As  such,  the  algorithm  proposed  is 
well-suited  for  real-time  DOT  image  reconstruction. 
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Abstract 

Diffuse  optical  tomography  (DOT)  poses  a  typical  ill-posed  inverse  problem 
with  a  limited  number  of  measurements  and  inherently  low  spatial  resolution. 
In  this  paper,  we  propose  a  hierarchical  Bayesian  approach  to  improve  spatial 
resolution  and  quantitative  accuracy  by  using  a  priori  information  provided  by 
a  secondary  high  resolution  anatomical  imaging  modality,  such  as  magnetic 
resonance  (MR)  or  x-ray.  In  such  a  dual  imaging  approach,  while  the 
correlation  between  optical  and  anatomical  images  may  be  high,  it  is  not 
perfect.  For  example,  a  tumour  may  be  present  in  the  optical  image,  but 
may  not  be  discemable  in  the  anatomical  image.  The  proposed  hierarchical 
Bayesian  approach  allows  incorporation  of  partial  a  priori  knowledge  about 
the  noise  and  unknown  optical  image  models,  thereby  capturing  the  function- 
anatomy  correlation  effectively.  We  present  a  computationally  efficient  iterative 
algorithm  to  simultaneously  estimate  the  optical  image  and  the  unknown 
a  priori  model  parameters.  Extensive  numerical  simulations  demonstrate 
that  the  proposed  method  avoids  undesirable  bias  towards  anatomical  prior 
information  and  leads  to  significantly  improved  spatial  resolution  and 
quantitative  accuracy. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Diffuse  optical  tomography  (DOT)  is  a  non-invasive  imaging  modality  that  makes  use  of 
the  light  in  the  near-infrared  (NIR)  spectrum  (Yodh  and  Chance  1995,  Hebden  et  al  1997, 
Arridge  and  Hebden  1997,  Intes  and  Chance  2005).  The  inverse  problem  in  DOT  involves 
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reconstruction  of  spatially  varying  absorption  and  scattering  properties  (O’Leary  1996,  Boas 
etal  2001,  Arridge  1999)  as  well  as  fluorophore  lifetime  and  yield  (Chang  etal  1997,  Hawrysz 
and  Sevick-Muraca  2000,  Eppstein  et  al  2002)  in  tissues  from  boundary  measurements. 
These  fundamental  quantities  can  be  utilized  to  obtain  tissue  oxy-  and  de-oxyhaemoglobin 
concentrations,  blood  oxygen  saturation,  water,  fat  and  contrast  agent  uptake  in  tissue  (Kincade 
2004).  The  unique  physiological  and  biochemical  information  offered  by  DOT  is  very  valuable 
for  practical  applications  such  as  breast  cancer  diagnosis  (Cerussi  et  al  2001,  Srinivasan  et  al 
2003,  Intes  et  al  2003),  cognitive  activity  monitoring  (Strangman  et  al  2002,  Villringer  and 
Chance  1997,  Chance  et  al  2003),  brain  tumour  and  haemorrhage  detection  (Hebden  et  al 
2004),  functional  muscle  imaging  (Quaresima  et  al  2003)  with  a  growing  list  of  applications 
in  fluorescence  tomographic  imaging  (Frangioni  2003,  Weissleder  and  Ntziachristos  2003). 

DOT  poses  a  typical  ill-posed  inverse  problem  with  a  large  number  of  unknowns 
and  a  relatively  limited  number  of  measurements.  This  necessitates  the  incorporation  of 
a  priori  information  into  the  inverse  problem  formulation  in  order  to  obtain  viable  solutions. 
Furthermore,  propagation  of  NIR  light  is  not  restricted  to  a  plane  owing  to  the  diffuse  nature 
of  photons  in  turbid  media,  which  results  in  poor  spatial  resolution.  To  tackle  the  ill-posed 
nature  of  the  inverse  problem  and  to  address  the  low  spatial  resolution  in  DOT,  a  number  of 
approaches  have  been  developed.  Bayesian  approach  has  been  suggested  to  incorporate  a  priori 
information  to  the  inverse  problem  formulation  (Oh  et  al  2002,  Milstein  et  al  2002,  Eppstein 
et  al  2002,  Guven  et  al  2002,  Ye  et  al  2001).  Introducing  penalty  functions  (Hielscher 
and  Bartel  2001)  and  uniform  (Paulsen  and  Jiang  1996,  Arridge  1993,  Jiang  et  al  1996, 
Yao  et  al  1997)  or  spatially  varying  regularization  terms  (Pogue  et  al  1999)  within  the 
regularization  framework  are  alternative  ways  to  incorporate  a  priori  information  into  the 
image  reconstruction  process.  In  all  these  studies,  no  other  imaging  modality  has  been  utilized 
to  infer  information  specific  to  the  medium  of  interest,  which  could  be  used  to  tune  the  prior 
information. 

1.1.  Related  literature 

Recently  several  research  groups  reported  development  of  hybrid  imaging  systems  combining 
optical  methods  with  high  resolution  anatomical  imaging  techniques.  These  include  a 
concurrent  x-ray  tomosynthesis-DOT  system  at  Massachusetts  General  Hospital  (Li  et  al 
2003),  MRI-DOT /DOS  (Diffuse  Optical  Spectroscopy)  systems  at  University  of  Pennsylvania 
(Intes  et  al  2002),  University  of  California  at  Irvine  (Gulsen  et  al  2003)  and  Dartmouth  College 
(Brooksby  et  al  2003)  and  ultrasound-DOT  system  at  University  of  Connecticut  (Zhu  et  al 
2003a).  These  multi-modality  developments  are  all  motivated  by  the  fact  that  DOT  offers 
unique  functional  information  (such  as  tissue  oxy-  and  deoxy-haemoglobin  concentrations) 
while  high  resolution  anatomical  imaging  modalities  provide  complementary  information  for 
disease  diagnosis  and  understanding  with  superior  localization  and  spatial  resolution.  Another 
incentive  comes  from  the  assumption  that  the  contrast  elements  provided  by  high  resolution 
imaging  modalities  correlate  well  with  the  optical  properties.  A  number  of  studies  lend 
support  to  this  assumption.  Ntziachristos  et  al  (2000)  have  reported  that  there  exists  a  good 
spatial  correlation  between  gadolinium  (Gd)-enhanced  MR  and  Indocyanine  Green  (ICG)- 
enhanced  DOT  images.  Cuccia  et  al  (2003)  have  also  shown  that  Gd-enhanced  and  Methyl 
Blue  (MB)-NIR  results  correlate  well  with  each  other  in  terms  of  perfusion  dynamics.  Merritt 
et  al  (2003a,  2003b)  presented  similar  observations  to  demonstrate  correlation  between  MR 
and  DOS  in  water  and  lipid  concentration  retrieval.  Furthermore  a  number  of  studies  have 
shown  that  incorporation  of  high  resolution  anatomical  images  as  a  priori  information  leads 
to  improved  diffuse  optical  image  reconstruction  (Dehghani  et  al  2002,  Ntziachristos  et  al 
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2002).  Pogue  and  Paulsen  (1998)  used  MR  images  to  generate  a  finite  element  mesh  to 
reconstruct  a  simulated  rat  cranium,  where  the  available  information,  about  the  structure  and 
optical  properties,  is  used  for  the  initial  guess  in  the  inversion  algorithm  similar  to  the  approach 
followed  by  Xu  et  al  (2002).  Schweiger  and  Arridge  (1999)  suggested  using  the  structural 
information  to  reconstruct  images  of  a  segmented  brain  model  at  a  low  resolution  level  in  order 
to  obtain  a  good  initial  guess  for  the  high-resolution  solution  of  the  same  problem.  Use  of  MR 
scans  has  been  employed  for  optical  breast  imaging  (Chang  et  al  1997,  Barbour  et  al  1995) 
where  the  ‘reference  medium’  was  obtained  from  accurate  optical  properties  of  the  tissue, 
with  the  anatomy  derived  from  MR  images.  Li  et  al  (2003)  reported  optical  breast  imaging 
results  guided  by  x-ray  mammography,  where  x-ray  contrast  was  assumed  to  be  proportional 
to  DOT  contrast.  X-ray  images  are  used  as  the  spatial  constraint  to  decompose  the  optical 
medium  into  two  major  sub-domains,  representing  the  region  of  interest  as  referred  to  the 
tumour  region  and  the  background,  respectively.  A  hybrid  L-curve  method  is  followed  for 
the  estimation  of  regularization  parameters  for  each  of  the  regularization  terms  corresponding 
to  the  sub-domains,  which  challenges  the  inverse  problem  computationally.  Brooksby  et  al 
(2003)  extended  the  idea  proposed  by  Schweiger  and  Arridge  (1999)  to  incorporate  the  initial 
low-resolution  optical  image  as  derived  from  MR  image  by  using  structural  information  and 
spatially  varying  regularization.  The  reported  results  are  encouraging;  however,  accurate 
quantification  of  the  tumour  region  is  possible  only  when  the  true  optical  heterogeneity  of 
tissue  distribution  is  included.  Therefore  in  this  approach,  the  overall  performance  relies  upon 
the  quality  of  the  initial  guess. 

1.2.  Proposed  method 

In  all  the  studies  referenced  above,  the  performance  of  the  DOT  image  reconstruction 
relies  on  the  assumption  that  the  correlation  between  the  anatomical  and  optical  images  is 
high.  However,  there  may  be  regions  in  the  optical  image  that  do  not  have  any  anatomical 
counterparts.  For  example,  a  tumour  may  be  apparent  in  the  optical  image,  but  may  not 
have  a  corresponding  signature  in  the  anatomical  image.  Furthermore,  average  optical 
coefficients  extracted  from  anatomical  images  may  be  significantly  different  from  the  true 
optical  coefficients  of  tissue.  As  a  result,  the  assumption  of  strong  optical-anatomy  correlation 
may  cause  undesirable,  erroneous  bias  in  optical  image  reconstruction.  Therefore,  more 
flexible  prior  models  are  needed  to  properly  represent  optical-anatomy  correlation.  For 
example,  when  the  average  optical  properties  extracted  from  anatomical  images  are  not 
reliable,  prior  image  model  should  provide  weaker  constraints  in  image  reconstruction.  The 
hierarchical  Bayesian  framework  affords  such  a  flexibility  in  designing  prior  image  and  noise 
models.  In  the  hierarchical  Bayesian  framework,  one  can  formulate  the  inverse  problem  in 
multiple  stages  where  each  stage  includes  information  about  the  unknown  parameters  of  the 
preceding  stage.  The  first  stage  of  the  hierarchy  includes  the  data  likelihood  and  the  first 
stage  of  the  image  prior,  which  comprise  statistical  models  for  the  noise  and  optical  image, 
respectively.  These  models  include  parameters  associated  with  noise  and  image  statistics, 
which  are  not  known  precisely  in  practice.  These  unknown  parameters  are  referred  to  as 
hyperparameters,  which  can  be  regarded  as  random  variables.  The  succeeding  stage  of 
the  hierarchical  formulation  incorporates  a  priori  information  about  the  hyperparameters  in 
the  form  of  prior  distributions — so  called  hyperpriors — defined  on  the  hyperparameters.  The 
incorporation  of  the  second  stage  concludes  the  design  of  the  two-level  hierarchical  noise  and 
image  models. 

In  this  paper,  we  consider  a  two-level  hierarchical  Bayesian  formulation  to  incorporate 
a  priori  anatomical  and  tissue  classification  information  into  the  DOT  image  reconstruction. 
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We  start  with  the  segmentation  of  the  high  resolution  image  and  classify  the  segmented 
image  into  sub-images  representing  major  tissue  types.  Based  on  the  tissue  label  information 
extracted  from  the  anatomical  image,  we  design  the  first  stage  of  the  prior  distribution  on  the 
unknown  optical  image  as  a  function  of  unknown  hyperparameters,  namely  the  image  mean  and 
standard  deviation.  Next,  we  design  the  data  likelihood  corresponding  to  the  parametric  noise 
model  with  an  associated  unknown  hyperparameter,  which  is  related  to  the  noise  variance.  The 
uncertainty  in  the  models  owing  to  the  unknown  hyperparameters  is  addressed  by  defining 
hyperpriors  on  the  hyperparameters,  which  constitutes  the  second  stage  of  the  hierarchical 
formulation.  The  hyperprior  on  the  noise-variance  related  hyperparameter  is  assumed  to  be 
uniform  so  as  to  not  constrain  its  value.  The  hyperpriors  on  the  hyperparameters  of  the  image 
model  are  formulated  with  the  aid  of  coregistered  tissue  classification.  Consequently,  the 
second  stage  of  the  image  prior  integrates  the  subjective  information  into  the  formulation, 
defining  the  extent  of  the  correlation  between  the  anatomical  and  optical  images.  As  a 
result,  a  priori  information  is  used  to  constrain  the  hyperparameters,  thereby  imposing  weaker 
constraints  on  the  optical  image.  We  refer  to  sections  2  and  3  for  a  detailed  discussion  of  the 
hierarchical  noise  and  image  models. 

Having  designed  the  hierarchical  noise  and  image  models,  we  formulate  the  joint 
distribution  of  the  measurements,  the  image  and  the  unknown  hyperparameters  associated  with 
the  noise  and  image  models.  In  order  to  estimate  the  hyperparameters,  we  adapt  the  linear 
conjugate  gradient  (CG)  algorithm  to  include  a  hyperparameter  estimation  step  followed  by  an 
image  update.  In  this  context,  we  apply  an  iterative  empirical  Bayesian  approach  to  estimate 
the  hyperparameters,  which  in  turn  gives  the  maximum  a  posteriori  (MAP)  estimates  of  the 
hyperparameters  at  each  CG  iteration  prior  to  the  image  update.  Hence,  the  noise  and  image 
models  are  accommodated  at  each  update  of  the  hyperparameters  along  with  the  solution 
process. 

We  perform  simulation  experiments  to  evaluate  the  performance  of  the  proposed 
hierarchical  Bayesian  formulation  and  hyperparameter  estimation  scheme.  Our  results  indicate 
that  hierarchical  Bayesian  approach  captures  the  function-anatomy  correlation  properly  and 
provides  improved  DOT  image  reconstruction  without  introducing  undesirable  bias  towards 
a  priori  anatomical  information.  Our  simulation  experiments  show  that  the  proposed  method 
provides  accurate  reconstruction  of  tumours  even  when  tumour  contrast  is  absent  in  the 
anatomical  image. 

1.3.  Organization  of  the  paper 

The  paper  is  organized  as  follows.  Section  2  defines  the  forward  model.  Section  3  provides 
background  on  the  hierarchical  Bayesian  formulation  of  the  inverse  problem  and  describes 
the  components  of  the  hierarchical  model.  Section  4  presents  the  iterative  algorithm  for  the 
simultaneous  estimation  of  the  optical  image  and  the  unknown  hyperparameters.  Section  5 
includes  numerical  experiments  to  validate  the  properties  of  the  proposed  approach.  Section  6 
summarizes  our  results  and  conclusion.  The  appendix  includes  the  derivation  of  the  estimation 
of  the  hyperparameters. 

2.  Forward  model 

In  the  NIR  region  of  the  electromagnetic  spectrum,  light  propagation  in  biological  tissue  can 
be  modelled  by  the  diffusion  approximation  to  the  radiative  transfer  equation.  The  diffusion 
equation  in  the  frequency  domain  is  given  by 

id) 

V  -D(r)  v  $(r)  -  pa(r)<p(r) - 0(r)  =  -A8(r  -  rf), 

c 
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where  <p(r)  represents  the  spatially  varying  optical  field  in  the  medium  £2  c  K3,  due  to  the 
point  source  AS(r  —rs)  located  at  r  =  rs.  co  denotes  the  modulation  frequency  of  the  source, 
c  is  the  speed  of  light  and  i  =  D(r)  is  the  spatially  varying  diffusion  coefficient  and 

pa  (r)  stands  for  the  spatially  varying  absorption  coefficient  of  the  medium. 

In  this  work,  we  focus  on  the  reconstruction  of  absorption  coefficients  of  the  medium. 
Therefore,  we  assume  that  the  diffusion  coefficient  of  the  medium  is  known  and  spatially 
invariant.  As  a  result,  the  following  diffusion  equation  given  in  frequency  domain  suffices  to 
define  the  forward  model: 

t  ico 

D  v2  <P(r)  -  pa(r)(j>(r) - cp(r)  =  -A8(r  -  rs).  (2) 

c 

We  have  employed  the  perturbation  approach  (O’Leary  et  al  1995,  Arridge  1995,  Kak 
and  Slaney  1988)  with  a  first-order  Rytov  approximation  to  solve  the  forward  problem  in  the 
frequency  domain  to  yield  a  system  of  linear  equations  after  the  discretization  of  the  medium 
£2  into  N  uniform  voxels  (Guven  et  al  2003a): 

y  =  Wx  +  £,  (3) 

where  y  is  the  measurement  vector,  W  is  the  Jacobian  based  on  the  Rytov  approximation 
computed  around  a  specified  homogeneous  background  f. ia(r)  =  pa o»  x  €  1R'V  denotes 
the  vector  of  differential  absorption  coefficients  Spa  of  the  medium  with  respect  to  the 
homogeneous  background  and  £  is  the  additive  noise  in  the  measurement  system.  Note 
that  recently  a  number  of  researches  (Ntziachristos  et  al  2002,  Li  et  al  2003,  Zhu  et  al  2003b  , 
Intes  et  al  2003)  have  reported  improved  DOT  reconstructions  for  clinical  images  based  on 
linearized  forward  model  using  high  resolution  anatomical  priors. 

3.  Hierarchical  Bayesian  formulation  of  the  inverse  problem 

We  approach  the  DOT  inverse  problem  from  a  Bayesian  perspective.  In  particular,  we  propose 
a  hierarchical  Bayesian  approach  to  effectively  capture  the  function-anatomy  correlation. 

We  formulate  the  posterior  distribution  of  the  unknown  image  and  compute  its  maximum 
a  posteriori  (MAP)  estimate  xMap;  that  is 

xMap  =  arg  max  {log  p  (x  [  y) } , 

X 

where  log  p(x|y)  is  the  posterior  distribution  of  the  unknown  image  x  given  the  measurements 
y.  Equivalently, 

£map  =  argmax{logp(y|x)  +  logp(x)}, 

X 

where  p{y\x)  is  the  data  likelihood  and  p(x)  is  the  prior  on  the  unknown  image  x. 

In  our  problem,  in  addition  to  the  optical  boundary  measurements  y,  we  also  have  the 
anatomical  tissue  label  information  C,  derived  from  the  a  priori  anatomical  image.  Therefore, 
the  MAP  estimate  can  be  modified  as 

&map  =  argmax{log  p(x\y,  C)} 

X 

=  arg  max  {log  p(y\x,  C)  +logp(x|C)},  (4) 

X 

where  p(y\x,C)  is  the  data  likelihood  and  p(x\C)  is  the  conditional  probability  density 
function  of  x  given  the  tissue  label  information  C. 

Given  the  forward  model  (3),  the  data  likelihood  is  governed  mainly  by  the  noise  statistics. 
Therefore  equation  (4)  reduces  to 

xMap  =  argmax{logp(y|x)  +logp(x|C)}. 


(5) 


2842 


M  Guven  el  al 


In  our  formulation,  the  noise  statistics  and  the  prior  distribution  are  governed  by  unknown 
model  parameters  A  and  Q,  respectively.  Here,  X  is  the  scalar  parameter  associated  with 
the  noise  variance  and  Q  is  the  vector  of  mean  and  variance  of  sub-images  that  correspond 
to  different  tissue  labels.  We  will  refer  to  these  parameters  as  hyperparameters.  In  order 
to  estimate  the  unknown  image  and  the  associated  hyperparameters,  we  consider  the  joint 
conditional  distribution  p(y,  x,  X,  Q\C): 

P(y ■  x,  X,  Q\C)  =  p{y ,  X\x)p(x,  Q\C),  (6) 

where  p{x,  Q\C)  is  the  conditional  hierarchical  prior  on  the  image  x  and  the  hyperparameters 
Q.  Equation  (6)  can  be  alternatively  expressed  as 

log  p{y,  x,  X,  Q\C)  =  log  p(y\x,  X)  +  log  p(X)  +  log  p(x\Q,  C)  +  log  p(Q\C),  (7) 

where  p{ X)  is  the  prior  distribution  on  X.  We  shall  refer  to  p{ X)  and  p(Q\C)  as  hyperpriors 
(Berger  1988).  In  this  representation,  p(x\Q,  C)  is  the  first-stage  prior  and  p(Q\C)  stands 
for  the  second-stage  prior.  In  the  following  sections,  we  will  discuss  how  the  data  likelihood 
and  the  hierarchical  prior  are  modelled. 


3.1.  The  data  likelihood  model 


The  measurement  vector  y  is  formed  as 

y  =  [y a  yn  •••  y{'D  ytl  •••  y«»  y(t  •••  y^f*  (g) 


where  S  is  the  number  of  sources,  D  is  the  number  of  detectors  and  F  is  the  number  of 
frequencies  associated  with  each  source.  The  total  number  of  measurements  is  then  equal  to 
P  =  S  x  D  x  F.  For  computational  efficiency,  we  limit  the  data  set  to  the  real  part  of  the 
measurements,  thus  y  eRp. 

Photon  detection  can  be  modelled  using  shot  noise  statistics,  which  originates  from 
Poisson  statistics.  With  a  sufficiently  large  number  of  detected  photons,  the  Poisson  statistics 
can  be  approximated  by  a  Gaussian  distribution,  with  a  variance  proportional  to  the  magnitude 
of  the  measurements  (Ye  et  al  2001,  Oh  et  al  2002).  Consequently,  we  model  the  data 
likelihood  in  equation  (6)  as 


p(y,k\x) 


K  |A?(A)|i/2 


exp 


(9) 


where  we  assume  a  non-informative  prior  for  A,  which  is  a  uniform  density  on  R1  (Berger 
1988).  In  the  above  distribution,  Af(A)  is  the  covariance  matrix  of  size  P  x  P,  K  is  the 
normalization  constant  and  ||z||^  =  zTAz.  Under  the  assumption  of  statistical  independence, 
Af  (A)  becomes  a  diagonal  matrix  of  the  form: 


Af(A)  =  AAj, 


'A  ol  0  0  •••  0 

0  A  o£  0  •  •  •  0 

0  0  0  0 

:  0 
0  0  0  0  Aor* 


(10) 


where  ct(2  is  equal  to  the  absolute  value  of  the  pth  measurement  and  the  unknown  parameter 
A  controls  the  scale  of  the  noise  covariance  matrix.  Therefore,  we  shall  refer  to  A  as  ‘noise 
scale’  for  the  rest  of  the  paper. 
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3.2.  Hierarchical  formulation  of  the  prior  distribution 


We  downsample  the  high  resolution  anatomical  image  to  match  its  resolution  with  that  of 
the  optical  image.  Next,  we  utilize  the  downsampled  anatomical  image  to  decompose  the 
optical  domain  into  M  non-overlapping  sub-images,  each  of  which  is  assumed  to  represent  a 
major  tissue  type.  For  instance,  for  breast,  these  tissues  types  can  be  parenchyma,  adipose 
and  tumour  (Ntziachristos  et  al  1999,  2000).  We  assume  a  Gaussian  distribution  for  each  sub¬ 
image,  with  unknown  mean  and  standard  deviation.  Thus,  the  first  stage  in  the  hierarchical 
prior  distribution  for  the  z  th  sub-image  is  given  by 


1  f  1  ,1 

p(Xi  | ^ ,  ad  =  - exp  -  —  2  Ha:,-  -  pi  ||2  , 

(2 no?)  ''  L  2 of  J 


i  =  1,2 . M 


(ID 


with  the  implicit  assumption  that  the  voxels  in  each  sub-image  are  statistically  independent. 
/x;  =  (f. ii  ■  •  ■  Pi)T  is  the  uniform  mean  value  vector  of  size  Nj  x  1,  where  Nj  stands  for 
the  number  of  voxels  in  the  ith  sub-image.  The  covariance  matrix  associated  with  the  ith 
sub-image  A*(er,)  =  crflNi-xN,  where  er,-  is  the  standard  deviation  of  each  voxel  in  the  ith 
sub-image  and  I *  n,  is  the  Nj  x  Nj  identity  matrix.  Assuming  that  the  sub-images  are 
statistically  independent,  the  first-stage  prior  of  the  image  given  the  tissue  label  information 
C  is 

p{x\Q,  C)  =  p{x \p,  a,  C ) 

1  ,  -  ,  . 
.  (12) 


'1  . 

' exp  --11*  -mIIajV)  - 


(2jr)w/2|Aj:(<T)|1/2 

where  Q  =  [p,  or],  p  is  the  vector  of  mean  values  assigned  to  the  sub-images  and  cr  is  the 
vector  of  standard  deviations  associated  with  the  sub-images,  that  is 
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(13) 


(14) 


and  A x(cr)  is  the  covariance  matrix  of  the  image  x 

f  N\xN\  0 


A*  (cr)  = 


0 


a}  I 


21N1xN1 


0 
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0 


0 


0  0 


(15) 


The  second-stage  prior  involves  incorporation  of  the  a  priori  information  into  the 
hierarchical  prior  distribution;  in  the  form  of  hyperpriors  defined  on  the  unknown 
hyperparameters  of  the  first  stage,  that  is  the  mean  and  standard  deviation  of  different  tissue 
types  in  the  optical  image.  Note  the  mean  values  of  different  tissue  types  are  specific  to  the 
unknown  optical  image  and  are  different  from  average  optical  values  of  tissues  available  in 
the  literature.  Nevertheless,  the  information  available  in  the  literature  can  be  used  to  design 
hyperpriors  on  the  unknown  mean  and  standard  deviation,  which  allows  effective  modelling 
of  the  uncertainty  in  the  prior  information. 

We  assume  a  Gaussian  distribution  for  the  mean  value  p  of  the  image: 

pWC)  =  exp  H'"  ‘  •  <16) 
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where  p  =  p\ and  pt  is  the  average  differential  absorption  of  the  i'th  tissue  type. 
AM(d)  =  A.t(<r)U=t?;  is  the  covariance  matrix  where  i?,-  stands  for  the  standard  deviation  of 
the  mean  value  of  each  voxel  in  the  ith  tissue  type,  for  i  =  1,2 , M. 

Similarly,  we  assume  a  Gaussian  distribution  for  standard  deviation  cr  of  the  image: 

pMC) = “p  B lk  ~  1  vw]  •  <17) 

where  pa  =  #*U=M«r-  a°d  A  a(y)  =  A^(<t)|  aj-Yi  for  i  =  1,2 , ,M.  Thus,  the  second-stage 
prior  p(Q\C)  becomes 

P(Q\C)  =  p(n,  <t\C)  =  p(p\C)p(<r\C) 

=  wwkg  “p  B*""  “  WV<«  +  lk  “  - 

(18) 

Having  designed  the  first-  and  second-stage  priors,  the  hierarchical  prior  distribution  in 
equation  (6)  becomes 

p(x,  QIC)  =  p(x,  p,  <r\C)  =  p(x \p,  a,  C)p(p\C)p(cr\C) 

1 

“  (27T)3^/2|AJC(0-)|1/2|A^(#)|l/2|Acr(y)|l/2 

xexP  +  +  •  (19) 

In  practice,  can  be  assigned  based  on  the  average  absorption  coefficients  of  tissue 
types  provided  in  the  literature  with  a  sufficiently  large  variance  $f.  Analogously,  the  mean 
value  of  a  can  be  extracted  from  the  error  bounds  of  the  average  optical  properties  of  tissue 
types  which  are  documented  in  the  literature.  See  for  example,  Mobley  and  Vo-Dinh  (2003) 
and  Cheong  et  al  (1990)  for  an  exhaustive  list  of  optical  coefficients  for  human  tissue  and 
fluids.  An  alternative  method  could  be  to  employ  a  low-resolution  estimate  of  the  optical 
image  (typically  a  least-squares  solution)  to  extract  approximate  values  for  these  parameters 
from  the  optical  data.  While  alternative  distribution  models  on  p  and  cr  can  be  considered,  we 
will  see  in  the  next  section  and  appendix  that  the  Gaussian  model  leads  to  a  computationally 
efficient  hyperparameter  estimation  scheme. 

4.  Image  reconstruction  and  hyperparameter  estimation 

Following  an  empirical  Bayesian  approach  (Berger  1988),  we  propose  an  iterative  algorithm  to 
estimate  both  the  optical  image  and  the  hyperparameters.  At  each  iteration,  the  MAP  estimates 
of  the  hyperparameters  are  computed  by  successively  maximizing  the  joint  distribution  with 
respect  to  each  hyperparameter.  The  hyperparameter  estimation  step  at  each  iteration  is 
followed  by  an  image  update. 

Substituting  equations  (9)  and  (19)  into  (6),  the  joint  probability  distribution  of  the 
measurements,  optical  image  and  the  hyperparameters  becomes 
p(y,  x,  X,  p,  cr\C)  =  p(y,  X\x)p{x\p,  cr,  C)p(p\C)p{cr\C) 

1 

“  (2tt )3^/2 1 Af  (A) |  !/2 1 A^  (or)  1 1/2 1 AM (r?)  1 1/2 1 ACT  (y ) I V2 

X  exp  -|(||2/  -  Wadl^+ll*  -  p\\\^(rT) 

+  \\p-M2A->m  +  \W-pJ2KHy))  • 


(20) 
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Let  '{'(x,  A,  p,  <t)  be  the  objective  functional  given  by 

^(x,  A.,  p,  cr)  =  -log  p{y,  x,  A,  p,  rr\C).  (21) 

Then,  the  image  reconstruction  becomes  an  optimization  problem  in  which  the  objective 
functional  '!'(*,  A,  p,  a)  has  to  be  minimized  with  respect  to  the  image  x  and  the 
hyperparameters  A,  p  and  cr. 

The  hyperparameter  estimation  problem  has  been  a  focus  of  both  the  statistical  and 
engineering  communities  and  many  procedures  have  been  suggested  to  date  (Mohammad- 
Djafari  1993, 1996,  Zhou  et  al  1994,  Molina  et  al  1999,  Utsugi  1997).  Since  the  optimization 
with  respect  to  the  optical  image  itself  is  a  computationally  intense  problem,  it  is  desirable  to 
keep  the  computation  complexity  of  hyperparameter  estimation  to  minimum.  Therefore,  we 
propose  an  iterative  algorithm  based  on  the  empirical  Bayesian  approach  that  successively 
estimates  the  hyperparameters.  The  hyperparameter  estimation  step  is  followed  by  the 
image  update  by  one  iteration  of  CG  algorithm,  applied  with  the  current  estimates  of  the 
hyperparameters. 

We  consider  an  alternating  minimization  scheme  for  the  estimation  of  the  mean  and 
standard  deviation,  where  the  current  estimate  of  one  of  the  parameters  is  used  to  estimate 
the  other.  This  approach  provides  a  computationally  efficient  solution  to  the  hyperparameter 
estimation  problem  with  only  O(N)  extra  operations  at  each  iteration,  where  N  is  the  number  of 
voxels.  The  estimation  of  the  noise  scale  A  is  independent  of  the  image  model  hyperparameters 
and  requires  O(N)  extra  operations  when  a  gradient-based  algorithm  (such  as  conjugate 
gradient  algorithm)  is  used.  As  a  result,  the  image  is  updated  based  on  the  current  estimates 
of  the  hyperparameters  with  O (N  P)  operations,  by  one  CG  iteration. 

A  pseudocode  describing  the  details  of  the  proposed  iterative  image  and  hyperparameter 
estimation  scheme  is  given  in  table  1.  The  outline  of  the  algorithm  is  as  follows:  the  image 
estimate  is  initialized  to  zero  vector  at  the  beginning  of  the  iterations.  At  each  iteration,  for 
the  given  image  update  x,  we  consider  x  as  the  hidden  variable  of  the  conditional  probability 
p(y\ A)  and  formulate  the  MAP  estimate  of  the  hyperparameter  A,  which  corresponds  to  the 
minimization  of  the  objective  functional  'I'(ai,  A,  p,  cr)  with  respect  to  A, 

Amap  4-  arg  min  ^(x ,  A ,  p,cr).  (22) 

Hyperparameters  of  the  image  model  are  estimated  in  a  similar  way.  We  formulate  the 
MAP  estimate  of  the  hyperparameter  p,  given  the  image  update  as  the  observations  and  the 
current  estimate  of  the  standard  deviation  as  the  hidden  variable  of  the  conditional  probability 
p(x\p,  C ).  This  corresponds  to  the  minimization  of  the  objective  functional  with  respect 
to  p, 

Amap  argmin  ^(x,  A, /z,  &).  (23) 

Similarly,  we  formulate  the  MAP  estimate  of  the  hyperparameter  cr,  given  the  image  update  as 
the  observations  and  the  current  estimate  of  the  mean  as  the  hidden  variable  of  the  conditional 
probability  p{x\cr,C).  This  is  equivalent  to  the  minimization  of  the  objective  functional  with 
respect  to  cr: 

o-map  4-  arg  min  'P(x,  A,  p,  a).  (24) 

a 

Note  that  the  estimate  A  is  attained  regardless  of  the  value  of  cr  and  p  and  vice  versa  (see  the 
appendix).  The  estimation  of  the  hyperparameters  is  followed  by  the  update  of  the  image  by 
one  CG  iteration 

a: MAP  4—  CGupdate{'h(x,  A,  p,  cr)}, 


L 


(25) 
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Table  1.  The  modified  conjugate  gradient  algorithm  with  the  proper  initializations  and  embedded 
hyperparameter  estimation  sub-routine. 

begin(Initialize) 

Image:  4®  =  0 

Hyperparameters:  X(0)  =  1;  £i  =  /t;  d/0)  =  t,-  for  i  =  1, 2, . . . ,  Mkj  >  0 
Gradient  vector:  </<°>  =  U-km 

Search  direction:  rf®  =  g<0) 

Termination  criterion:  e 
Iteration  counter:  n  =  0 
end(Initialize) 
repeat 

begin{Update  Image) 

Exact  line  search:  =  arg  mina>o  6^) 

4(b+,)  =  £(b>  +a(n)Sn) 

end{Update  Image)  J 

beginfEstimate  Hyperparameters) 

X<"+l>  <—  arg  minx  4<  (x^"\  X,  a^n>) 

*)  <—  argmin^  X^n\  /x,  a(-n)) 

d<"+1>  -r-  arg minCT  *(*<">,  A<">,  A(n+1),  o') 
end(Estimate  Hypeiparameters) 
begin(Update  Search  Direction) 
p(n+l)  _  -£2(4(n+1>) 

^"l)=max(^yy>),0) 

rf(n+D  _  fl(n+l)  +  p(n+])d(n) 
end(Update  Search  Direction) 
n  =  n  +  1 

until  ||q'('i)S2(£(',))||  <  £ 


where  the  step  length  is  computed  by  the  exact  line  search  (Nash  and  Sofer  1996)  and  Polak- 
Ribiere-Polyak  method  (Polak  andRibiere  1969,  Polyak  1969)  is  used  to  calculate  fi  parameter 
of  the  CG  method.  The  derivation  of  the  hyperparameters  is  given  in  the  appendix  in  detail. 
The  computational  complexity  of  the  proposed  image  reconstruction  and  hyperparameter 
estimation  algorithm  is  shown  in  table  2. 

Note  that  the  proposed  simultaneous  image  reconstruction  and  hyperparameter  estimation 
algorithm  can  be  viewed  as  a  variation  of  the  alternating  minimization  algorithm  (Csiszar  and 
Tusnady  1984),  where  the  minimization  of  the  objective  functional  W(x,  A,  /x,  a)  with  respect 
to  the  image  is  replaced  by  one  CG  iteration  that  leads  to  the  update  of  the  image.  Similar 
approaches  can  be  found  in  Mohammad-Djafari  (1993, 1996)  and  Milstein  et  al  (2002).  The 
empirical  approach  proposed  in  this  work  is  asymptotically  efficient  and  comparable  with  the 
hierarchical  analysis  (Berger  1988,  Molina  et  al  1999)  provided  that  the  number  of  observations 
(Ni  for  the  sub-images,  P  for  the  measurements)  is  large  (Berger  1988).  Alternative  approaches 
for  hyperparameter  estimation  include  ML-Type  II  (Berger  1988)  and  marginalized  ML  type 
estimation  schemes,  which  do  not  incorporate  hyperpriors.  Extension  of  these  approaches 
to  marginalized  MAP  estimation  requires  integration  over  the  multi-dimensional  image  x 
(Molina  et  al  1999,  Galatsanos  et  al  2002),  which  may  result  in  increased  computational 
complexity. 

5.  Results 

We  perform  three  sets  of  experiments  to  evaluate  the  performance  of  the  proposed  method.  For 
each  of  the  experiments,  we  used  a  finite  difference  code  to  simulate  the  optical  measurements. 
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Table  2.  Computation  complexity  of  the  algorithm  described  in  table  1 .  N  is  the  number  of  voxels 
and  P  is  the  number  of  measurements.  Since  the  noise  scale  estimate  X  is  used  in  the  calculation 
of  a,  this  saves  O(iVP)  number  of  multiplications  and  O (N P)  number  of  additions.  Thus, 
the  estimation  sub-routine  does  not  affect  the  overall  computational  complexity  of  the  original 
conjugate  gradient  algorithm. 


Operation 

Number  of 
multiplications 

Number  of 

additions 

a  computation 
(exact  line  search) 

0(N  P) 

O  (NP) 

ft  parameter  calculation 

O(N) 

O  (N) 

Estimation  of  X 

O  (NP) 

O  (NP) 

Estimation  of  p 

O(N) 

O  (A0 

Estimation  of  <r 

O(N) 

O(N) 

Gradient  calculation 
and  image  update 

O  (NP) 

O  (NP) 

Total 

O  (NP) 

O  (NP) 

3%  of  the  mean  value  of  measurements  was  added  to  the  measurement  vector  y  in  each 
experiment. 

In  the  first  set  of  experiments,  the  objective  is  to  evaluate  whether  the  optical  image 
reconstruction  is  biased  towards  the  average  optical  coefficients,  which  are  used  to  formulate 
the  hyperprior  defined  on  the  mean  value.  Second  set  of  experiments  demonstrate  how  well 
the  hierarchical  Bayesian  formulation  captures  the  correlation  between  the  optical  and  the 
a  priori  anatomical  image.  In  the  last  set  of  experiments,  we  evaluate  the  proposed  method 
using  optical  data  simulated  from  an  MR  breast  image.  We  show  that  the  a  priori  information 
improves  the  image  reconstruction  and  does  not  lead  to  an  erroneous  bias  towards  the  a  priori 
information. 


5.1.  Simulation  experiment  I 

A  priori  selection  of  parameters  in  the  assumed  image  and/or  noise  models  may  bias  the  optical 
image  reconstruction  in  an  undesirable  way.  The  hierarchical  Bayesian  formulation  and  the 
empirical  hyperparameter  estimation  scheme  proposed  in  this  paper  avoids  such  undesirable 
results  by  incorporating  dynamic  image  and  noise  models  in  the  problem  formulation.  In 
this  experiment,  we  show  that  the  proposed  hyperparameter  estimation  approach  is  relatively 
insensitive  to  average  optical  values  used  to  design  the  hyperprior  p{p\C),  defined  on  the 
mean  value  of  the  image. 

We  consider  a  square  heterogeneity  with  a  mean  absorption  coefficient  of  0.071  cm-1 
embedded  in  a  background  with  a  mean  absorption  coefficient  of  0.04  cm-1  as  shown 
in  figure  1.  The  diffusion  coefficient  of  both  the  heterogeneity  and  the  background  is 
assumed  equal  and  set  to  D  =  0.033  cm.  We  consider  a  transmission  geometry  and 
distribute  19  sources  and  19  detectors  on  opposite  sides  to  yield  a  total  of  722  measurements, 
collected  at  two  frequencies,  that  is  100  and  200  MHz.  We  evaluate  the  Jacobian  at  pao  = 
0.04  cm"1  and  perform  200  experiments  with  the  proposed  hierarchical  Bayesian  formulation 
and  hyperparameter  estimation  scheme.  At  each  experiment,  p(p\C)  is  formulated  such  that 
Pi  +  pa0  for  the  square  inclusion  and  for  the  background  are  drawn  randomly  from  a  uniform 
distribution  with  lower  and  upper  bounds  (0.038,  0.114)  and  (0.02,  0.06),  respectively  (Note 
that  pi  value  is  used  in  the  formulation  of  the  objective  functional).  The  associated  standard 
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Figure  1.  The  medium  used  to  simulate  the  optical  data  for  the  first  experiment.  A  square  absorber 
is  embedded  in  an  almost  homogeneous  background. 
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Figure  2.  The  average  absorption  values  ( p,a )  in  the  reconstructed  sub-images  for  randomly  drawn 
average  absorption  values  /x„.  Results  are  shown  for  200  trials.  Note  that  the  (ia  values  estimated 
for  the  square  inclusion  (an  average  value  of  0.068)  are  very  close  to  the  actual  value  of  0.071. 
(a)  The  average  value  in  the  reconstructed  square  heterogeneity  versus  the  average  absorption 
value  drawn  for  the  square  heterogeneity  is  shown,  (b)  The  average  value  [la  in  the  reconstructed 
background  sub-image  versus  the  average  absorption  value  p,a  drawn  for  the  background  is  shown. 


deviation  is  sufficiently  large  and  set  to  t>;  =  6 (//.,•  +  /x„o).  Similarly  p(cr\C)  is  formulated 
such  that  =  0.4(/x,  +  /to0)  and  =  15/tff,. 

This  simulation  study  demonstrates  that  the  proposed  method  provides  effective  means 
to  constrain  image  reconstruction  without  biasing  the  solution.  Figure  2  shows  the  average  of 
the  estimated  absorption  values  /t a  (i.e.  sr,  +  p,a 0)  versus  the  assigned  hyperparameter  values 
jla  (i.e.  jXi  +  /xao)  for  the  square  inclusion  and  the  background,  respectively.  We  observe  that 
quantitative  accuracy  is  achieved  even  in  the  extreme  cases  and  the  reconstruction  for  the 
background  is  almost  insensitive  to  the  assigned  hyperparameters. 

5.2.  Simulation  experiment  II 

In  this  experiment,  we  examine  a  case  where  the  heterogeneity  is  present  in  the  a  priori 
anatomical  image  but  not  in  the  optical  image  (see  figure  3).  Based  on  the  anatomical 
template,  the  optical  image  was  segmented  into  two  sub-images,  one  corresponding  to  the 
background  and  the  other  corresponding  to  the  two  inclusions,  which  were  assumed  to  have 
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Figure  3.  The  optical  image  (a)  and  the  anatomical  counterpart  (b).  Note  the  additional  absorber 
(square  inclusion)  indicated  by  the  anatomical  prior  which  does  not  exist  in  the  optical  map. 

Table  3.  The  parameter  set  used  in  the  inverse  problem  formulations  for  the  simulation  experiment 
II,  displayed  for  the  two  sub-images:  the  inclusions  and  the  background.  pao  =  0.04  and 
D  =  0.33  cm  for  this  experiment.  This  experiment  was  performed  with  the  same  source-detector 
configuration  as  in  experiment  1. 


Hyperprior  parameters 

(fii  +  Pad,  #i) 

(Pa,,Yi) 

Both  inclusions  in 

(0.08, 0.04) 

(0.03, 0.45) 

the  anatomical  image 
Background 

(0.04, 0.02) 

(0.02, 0.30) 

the  same  average  optical  coefficients.  For  comparison,  we  also  considered  the  maximum 
likelihood  (ML)  approach  for  the  inverse  problem  formulation.  The  ML  approach  estimates 
the  optical  image  based  on  the  data  likelihood  model  given  in  section  3.1.  This  formulation 
does  not  incorporate  any  a  priori  information  about  the  image.  Nevertheless,  the  noise  scale 
A.  is  unknown  and  has  to  be  estimated  as  described  in  section  4.  Note  that  the  ML  (no  prior) 
formulation  is  in  fact  regularized  by  the  stopping  criterion  of  the  conjugate  gradient  algorithm 
used  in  the  minimization  of  the  resulting  objective  functional.  We  list  the  parameter  set  used 
in  the  hierarchical  prior  design  in  table  3. 

The  reconstructed  images  are  shown  in  figure  4.  The  true  mean  values  and  the  sample 
average  of  the  estimated  absorption  coefficients  of  the  heterogeneity  on  the  left  and  on  the  right 
and  of  the  background  are  shown  in  table  4.  Even  though  the  anatomical  image  indicates  a 
heterogeneity  on  the  right,  the  hierarchical  Bayesian  formulation  leads  to  a  qualitatively  good 
reconstruction.  The  ML  estimate  of  the  image  detects  the  rectangular  inclusion,  but  suffers 
from  low  resolution  and  lacks  accuracy  in  the  reconstructed  value  of  the  absorption  coefficient 
of  the  rectangular  absorber. 

5.3.  2D  experiment  with  MR-simulated  data 

We  used  the  T1 -weighted  MR  breast  image  from  Ntziachristos  etal(  1999)  to  design  a  realistic 
optical  breast  model  (figure  5).  The  MR  breast  image  was  segmented  into  parenchyma  and 
adipose  layers  by  applying  a  simple  thresholding  algorithm  with  respect  to  the  MR  image 
intensity  values.  Next,  a  tumour  corresponding  to  an  infiltrating  ductal  carcinoma  revealed 
by  Gd-DTPA  (gadolinium-diethylenetriamine  pentaacetic  acid)  enhancement  was  inserted 
(shown  in  figure  5  as  well).  Each  sub-region  was  assigned  an  absorption  value  as  indicated 
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(a)  (b) 


Figure  4.  The  reconstructed  images  as  a  result  of  hierarchical  Bayesian  formulation  (a),  and 
ML  (no  prior)  formulation  (b).  ML  solution  is  quantitatively  inaccurate  and  suffers  from  low 
spatial  resolution.  Hierarchical  Bayesian  formulation  leads  to  an  image  estimate  with  quantitative 
accuracy. 


Figure  5.  The  original  MR  breast  image  with  an  artificial  tumour  inserted. 


Table  4.  The  actual  mean  values  and  the  mean  of  the  reconstructed  sub-images  as  a  result  of 
maximum  likelihood  (no  prior)  and  hierarchical  Bayesian  formulations.  Sub-images  are  defined 
on  the  sub-domains  as  indicated  by  the  anatomical  image  shown  in  figure  3(b). 


Sub-domain 

True  mean 

value 

Maximum 

likelihood 

Hierarchical 

Bayesian 

Rectangular  inclusion 

0.076 

0.061 

0.071 

Square  inclusion 

0.040 

0.037 

0.044 

Background 

0.040 

0.041 

0.040 

in  Ntziachristos  et  al  (1999)  (^dipost'  =  0.03  cm-1,  /t4Parench>'ma  =  0.06  cm-1,  ^umour  = 
0.09  cm-1)  to  obtain  an  initial  template  (figure  6(b)).  To  simulate  a  corresponding  optical 
image,  zero  mean  Gaussian  noise  was  added  prior  to  filtering  the  image  by  a  low-pass  filter. 
The  resulting  optical  image  is  shown  in  figure  6(b).  Note  the  quantitative  and  spatial  mismatch 
along  the  boundaries  and  especially  within  the  tumour.  The  homogeneous  diffusion  coefficient 
of  the  medium  was  set  to  0.042  cm.  Nine  frequencies  ranging  from  0  to  244  MHz  were  used 
to  obtain  729  measurements  with  nine  sources  and  nine  detectors  positioned  along  the  x-axis 
on  opposite  sides.  The  optical  medium  was  uniformly  discretized  into  90  pixels  along  the 
x-axis  and  60  pixels  along  the  y-axis  leading  to  a  total  of  5400  lxl  cm2  pixels. 
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(a)  (b) 


Figure  6.  The  anatomical  and  optical  images  are  shown  on  the  left  and  right,  respectively.  Note 
the  spatial  mismatch  between  the  two  images. 


(c)  (d) 


Figure  7.  The  hierarchical  Bayesian  reconstruction  of  the  optical  image  (a)  using  the  anatomical 
template  shown  in  figure  6(a)  for  the  design  of  the  hierarchical  image  model.  Part  (c)  shows  the 
image  that  zooms  into  the  tumour  region  in  the  optical  image  shown  in  (a).  The  ML  estimate  of  the 
entire  image  and  the  sub-image  focusing  the  tumour  region  are  shown  in  (b)  and  (d),  respectively. 
The  rectangular  box  in  the  figures  shows  the  actual  location  of  the  tumour. 


We  performed  two  types  of  experiments  to  test  the  performance  of  the  proposed 
hierarchical  Bayesian  approach  for  this  problem. 

(i)  Tumour  present  both  anatomically  and  optically.  In  this  experiment,  the  template 
extracted  from  the  anatomical  image  shown  in  figure  6(a)  was  used  to  design  the  hierarchical 
image  prior.  As  a  result,  the  optical  image  was  segmented  into  three  sub-images  each  of 
which  corresponded  to  the  labelled  images  in  the  anatomical  image  as  shown  in  figure  6(a). 
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Figure  8.  The  mean  value  and  the  standard  deviation  estimates  for  each  sub-image  (sub-images 
are  determined  by  the  anatomical  image)  and  noise  scale  estimate  versus  iteration  number  are 
shown.  The  thick  solid  line  shows  the  estimated  values.  The  constant  solid  line  in  (a)  shows  the 
actual  mean  value  and  the  constant  dashed  line  in  (a)  shows  the  assigned  mean  value  (/x,-  +  ixa o) 
used  in  the  design  of  the  hyperprior  defined  on  the  sub-image  means.  The  actual  mean  values  in 
these  sub-regions  are  0.032,  0.058  and  0.076,  respectively,  (a)  The  estimated  mean  values  /x,-  for 
each  sub-image  versus  iteration  number  are  shown.  The  estimates  for  the  parenchyma,  adipose 
and  tumour  sub-images  are  given  at  the  top,  middle  and  bottom,  respectively,  (b)  The  estimated 
standard  deviation  rf,  for  each  sub-image  versus  iteration  number  are  shown  on  the  left.  The 
estimates  for  the  parenchyma,  adipose  and  tumour  sub-images  are  given  at  the  top,  middle  and 
bottom,  respectively.  The  noise  scale  estimate  A  versus  iteration  number  is  shown  on  the  right. 


In  the  design  of  the  hyperprior  defined  on  the  mean  (i.e.  p(p,\C)),  values  that  are  significantly 
different  from  the  actual  mean  of  the  sub-images  were  used.  Thus,  this  experiment  evaluates 
the  robustness  of  the  proposed  method  when  the  true  statistics  of  the  optical  image  are 
significantly  different  from  the  statistics  extracted  from  the  prior  anatomical  image. 

The  reconstructed  image  and  the  sub-image  zoomed  into  the  tumour  region  are  shown  in 
figures  7(a)  and  (c),  respectively.  For  comparison,  the  ML  estimate  of  the  image  is  shown  in 
figures  7(b)  and  (d).  The  simulation  results  show  that  hierarchical  Bayesian  approach  leads  to 
qualitatively  better  results  and  resolves  the  tumour  more  accurately. 

In  figure  8,  the  estimates  of  the  hyperparameters  associated  with  the  noise  and  image 
models  are  given  as  a  function  of  the  iteration  number.  Note  that  the  mean  value  estimates 
for  each  sub-image  converge  to  actual  values,  even  though  the  corresponding  assigned 
hyperparameters  regarding  the  mean  value  deviate  from  the  true  average  optical  values  by 
at  least  15%  (see  table  5).  The  experiment  also  demonstrates  that  the  initialization  of  the 
hyperparameters  does  not  have  any  effect  on  the  performance  of  the  estimation  (figure  8). 

(ii)  Tumour  present  optically  but  not  anatomically.  In  this  experiment,  we  removed  the 
tumour  region  from  the  template  extracted  from  the  prior  anatomical  image,  but  kept  it  in  the 
optical  image  as  shown  in  figures  9(a)  and  (b).  As  a  result,  the  optical  image  was  segmented 
into  two  sub-images.  The  objective  of  this  experiment  is  to  evaluate  how  well  the  proposed 
method  reconstructs  optical  tumours  when  they  are  not  anatomically  present. 

The  reconstructed  images  for  this  experiment  are  given  in  figures  10(a)  and  (c), 
respectively.  The  ML  estimate  of  the  image  is  given  in  figures  10(b)  and  (d).  We  observe 
that,  even  though  there  is  a  significant  mismatch  between  the  optical  image  and  the  anatomical 
counterpart  in  the  tumour  region,  the  hierarchical  Bayesian  formulation  leads  to  a  qualitatively 
better  reconstruction  than  the  ML  approach,  even  around  the  tumour.  Furthermore,  the  tumour 
is  better  localized  as  compared  to  the  ML  solution  and  is  not  biased  towards  the  a  priori 
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(a)  (b) 

Figure  9.  The  anatomical  template  (a)  and  the  original  optical  image  (b).  Note  that  the  tumour  is 
not  anatomically  present. 


Figure  10.  The  hierarchical  Bayesian  reconstruction  of  the  optical  image  (a)  using  the  anatomical 
template  shown  in  figure  9(a)  for  the  design  of  the  hierarchical  image  model.  Part  (c)  shows  the 
image  that  zooms  into  the  tumour  region  in  the  optical  image  shown  in  (a).  The  ML  estimate  of  the 
entire  image  and  the  sub-image  focusing  the  tumour  region  are  shown  in  (b)  and  (d),  respectively. 
The  rectangular  box  in  the  figures  shows  the  actual  location  of  the  tumour. 

anatomical  image.  The  error  in  the  localization  of  the  tumour  can  be  attributed  to  the  source- 
detector  geometry.  The  propagation  of  light  along  the  y -direction  results  in  a  smoothing  effect 
on  the  optical  image  along  the  y-direction.  This  effect  is  enhanced  near  source  and  detectors 
due  to  the  behaviour  of  the  solution  of  the  diffusion  equation.  The  vertically  smoothing  effect 
can  be  observed  in  the  ML  estimate  of  the  image  more  apparently  (figures  10(b)  and  (d)).  The 
smoothing  effect  can  be  suppressed  by  incorporation  of  a  priori  information  for  the  tumour 
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(b) 


Figure  11.  The  mean  values  of  the  reconstructed  sub-images  versus  iteration  number  (a).  The 
estimated  values  of  the  standard  deviation  of  the  two  sub-images  and  the  noise  scale  X  are  shown  in 
(b).  (a)  The  average  value  ( p,a  in  each  reconstructed  sub-image  versus  iteration  number  are  shown. 
The  estimates  for  the  parenchyma,  adipose  and  tumour  sub-images  are  given  at  the  top,  middle 
and  bottom,  respectively.  The  domains  of  the  sub-images  correspond  to  the  domains  in  the  optical 
image  shown  in  figure  9(a).  (b)  The  estimated  standard  deviation  dj  for  each  sub-image  versus 
iteration  number  are  shown  on  the  left.  The  estimates  for  the  parenchyma,  adipose  sub-images  are 
given  at  the  top  and  bottom,  respectively.  The  sub-images  correspond  to  the  sub-images  defined  by 
the  anatomical  template  shown  in  figure  9(b).  The  noise  scale  estimate  X  versus  iteration  number 
is  shown  on  the  right. 

Table  5.  The  actual  mean  of  the  absorption  values  in  each  sub-image  and  the  parameter  set  used 
in  the  inverse  problem  formulations  for  the  MR-simulated  experiments  I  and  II.  /r„o  =  0.0439  for 
this  experiment.  N/A  stands  for  ‘not  applicable’. 


Sub-images 

Parenchyma 

Adipose 

Tumour 

The  first 

(£;  +/A.o,tfi): 

(0.038, 0.23) 

(0.05,  0.3) 

(0.09, 0.54) 

experiment 

.  Xi): 

(0.015, 0.228) 

(0.02,  0.3) 

(0.036,  0.54) 

The  second 

(Ai  +  MaO.  t?i): 

(0.03, 0.18) 

(0.06,  0.36) 

N/A 

experiment 

(Mcr, ,  Vi): 

(0.012,0.18) 

(0.024, 0.36) 

N/A 

(£aCtUal  +  Pat))'- 

0.032 

0.058 

0.076 

as  in  case  (i),  where  the  tumour  is  better  resolved  (figure  7(c)).  Further  improvement  can 
be  achieved  by  employing  sources  and  detectors  positioned  along  the  y-axis  as  well  as  along 
x-axis. 

In  figure  1 1 ,  the  average  absorption  value  of  each  reconstructed  sub-image  versus  iteration 
number  is  given.  The  sub-images  correspond  to  those  as  indicated  by  the  actual  optical  image 
shown  in  figure  6.  Note  that  the  mean  value  of  the  reconstructed  image  in  the  tumour  region 
converges  to  the  actual  value  even  though  the  anatomical  image  asserts  that  no  tumour  exists. 

The  set  of  parameters  used  in  the  design  of  hyperpriors  for  these  experiments  and  the 
actual  mean  of  absorption  values  for  each  sub-image  are  shown  in  table  5. 

6.  Conclusion 

In  this  work,  we  formulated  the  inverse  DOT  problem  within  a  hierarchical  Bayesian 
framework  where  the  hierarchical  prior  distribution  is  based  on  the  a  priori  information 
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extracted  from  a  secondary  high-resolution  anatomical  image.  Instead  of  directly  constraining 
the  optical  image  with  the  anatomical  prior,  we  incorporated  the  a  priori  information  in  the 
form  of  hyperpriors  to  impose  constraints  on  the  unknown  hyperparameters  of  the  image 
and  noise  models.  We  proposed  a  computationally  efficient  iterative  algorithm,  based  on  an 
empirical  Bayesian  approach,  to  simultaneously  estimate  the  optical  image  and  the  unknown 
hyperparameters.  We  tested  the  proposed  approach  in  three  different  simulation  experiments. 
Numerical  experiments  demonstrate  that  the  proposed  approach  improves  the  spatial  resolution 
and  quantitative  accuracy  of  optical  images.  Our  study  shows  that  the  hierarchical  Bayesian 
approach  provides  an  effective  framework  to  capture  the  correlation  between  optical  and 
anatomical  images. 

The  proposed  hierarchical  Bayesian  formulation  can  be  extended  to  incorporate  spectral 
a  priori  information  (Intes  et  al  2004).  Finally,  we  note  that  the  results  are  based  on 
the  linearized  forward  model  around  a  homogeneous  background.  However,  the  proposed 
hierarchical  Bayesian  formulation  and  the  iterative  optical  image  and  hyperparameter 
estimation  scheme  can  be  adapted  to  the  nonlinear  inverse  DOT  problem,  wherein  the  Jacobian 
of  the  forward  model  is  iteratively  updated. 
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Appendix 


The  minimization  of  the  objective  function  '!'(*,  A,  p,  cr)  with  respect  to  the  noise  scale  A 
given  the  updated  image  estimate  x  results  in  the  ML  estimate  A  of  the  noise  scale: 

X=^\\y-W<t\\2A.i.  (A.l) 

F  y 

Note  that  the  estimate  A  is  independent  of  the  hyperparameters  associated  with  the  image 
model.  On  the  other  hand,  the  minimization  of  the  objective  functional  with  respect  to  p 
yields  the  MAP  estimate  of  the  mean.  In  order  to  find  an  estimate  for  the  vector  p,  we  make 
use  of  the  probability  density  function  formulation  for  each  individual  sub-image  and  rewrite 
the  objective  function  with  Pi  dependent  terms,  given  the  sub-image  estimate  4,  and  the 
current  estimate  of  the  standard  deviation  b, : 

(**/)  =  ^  "  A/  II2  +  ^2  II /*«  -  A,  II2-  (A.2) 

Minimization  with  respect  to  yields  the  estimate  p,  for  the  mean  value  of  the  ith  sub-image, 

£»,) 

V,/lA ,  (A.3) 


__L| 

(  v?  A  e? 

1  t  Y ,  i  -4-  -  * 

Ni  ' 

_ L  j 

(  y , 

Ni  1 

where  x,*  denotes  the  kth  voxel  in  the  ith  sub-image  and  p =  jli  stands  for  the  assigned 
mean  value  of  the  image  mean  for  the  Ath  voxel  in  the  ith  sub-image. 


•4 


2856 


M  Guven  et  al 


Following  a  similar  procedure,  we  make  use  of  the  hierarchical  prior  formulation  for  each 
individual  sub-image  and  rewrite  the  objective  functional  with  ct,- -dependent  terms,  given  the 
sub-image  estimate  and  the  current  estimate  of  the  mean  A, : 

arg min  (or,)  =  log  (of')  +  ^ p,-  -  A,  II2  +  ^ Ik,  -  Ma,  II2-  (A.4) 

where  rr,  =  (o,-  •  •  •  o,)r.  After  taking  the  derivative  of  the  above  expression,  the  estimate  for 
the  sub-image  standard  deviation  a,-  satisfies  the  following  equation: 

Ni&l  -  Pi  -  A ,  II2  +  -  Mff,)  =  0.  (A.5) 

yf 

This  is  a  fourth-order  equation  in  d;,  which  is  rather  difficult  to  solve.  In  order  to  simplify  the 
solution,  we  make  use  of  the  following  approximation: 

ixaidf  ^  pi',  (A.6) 


and  equation  (A.5)  becomes  a  quadratic  equation  of  the  estimate  &j 

+  Ni&f  -  ^pi  -  A,  II2  +  Ni  =  0. 

Then  the  positive  value  for  the  variance  estimate  &f  is  equal  to 

-  -Nt  +  VA  A 

6?  = - L  >  0, 

'  2  Nt/y? 

where  the  discriminant  A  in  the  equation  can  be  evaluated  as 


=  Nf+4^  ^P,  -  A,  II2  +  N,  >  0. 


(A.7) 


(A.8) 


(A.9) 


Note  that  in  the  limiting  case,  for  yf  oo,  the  sub-image  variance  estimate  af  converges  to 


lim  <7, 2  =  lim 


-Ni+y/E  Pi  “Ail 


i  =  1,2,  ...,Af, 


jf4“oo  2  Ni/yf  Ni  ' 

which  corresponds  to  the  ML  estimate  of  the  variance  (Guven  et  al  2003b). 


(A.  10) 
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Abstract 

Diffuse  optical  tomography  is  a  typical  inverse  problem  plagued  by 
ill-condition.  To  overcome  this  drawback,  regularization  or  constraining 
techniques  are  incorporated  in  the  inverse  formulation.  In  this  work,  we 
investigate  the  enhancement  in  recovering  functional  parameters  by  using 
physiological  and  spatial  a  priori  constraints.  More  accurate  recovery  of 
the  two  main  functional  parameters  that  are  the  blood  volume  and  the  relative 
saturation  is  demonstrated  through  simulations  by  using  our  method  compared 
to  actual  techniques. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Probing  human  tissue  with  near  infrared  (NIR)  light  is  emerging  as  a  new,  promising  imaging 
modality.  The  strength  of  this  new  technology  relies  on  its  ability  to  reveal  the  functional  state 
of  deep  tissue.  To  date  the  main  applications  of  NIR  technologies  are  functional  brain 
imaging  (Strangman  et  al  2002,  Villringer  and  Chance  1997),  muscle  imaging  (Lin  etal 
2002)  and  optical  mammography  (Colak  etal  1999,  Franceschini  etal  1997,  Grosenick  etal 
2003,  Intes  et  al  2003,  Jiang  et  al  2002,  McBride  et  al  2001,  Tromberg  et  al  2000).  The  last 
application  is  known  to  be  of  growing  interest  in  both  the  research  and  medical  community 
(Cerussi  and  Tromberg  2003). 

Optical  mammography  aims  to  retrieve  important,  local  physiological  parameters  that 
are  the  blood  content  and  the  relative  oxygenation  of  the  blood.  These  two  functional 
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parameters  are  known  to  provide  a  means  of  discriminating  between  healthy  and  diseased 
tissues.  The  blood  content  relates  to  the  angiogenesis  level  of  the  tumourous  mass,  and  the 
relative  oxygenation  to  its  hypermetabolic  state.  These  two  functional  signatures  correlate  to 
malignancy. 

Pre-clinical  data  strengthen  these  assumptions.  However,  difficulties  arise  due  to  the 
specific  nature  of  the  light  propagation.  In  this  spectral  window,  the  light  is  strongly 
diffused  leading  to  a  relatively  poor  resolution  for  thick  tissue  investigation.  Also,  the 
breast  is  by  nature  a  heterogeneous  organ  with  a  complex  spatial  distribution  of  optical 
relevant  chromophores.  In  the  NIR  spectral  range,  four  chromophores  are  accountable  for 
the  absorption  (Cerussi  et  al  2001).  These  chromophores  are  the  oxy-  ([Hb02]),  deoxy- 
haemoglobin  ([Hb]),  the  water  ([H20])  and  the  lipids  ([Li]).  The  first  two  chromophores  are 
providing  an  insight  into  the  tissue  functional  state.  The  last  two  chromophores  are  linked  to 
the  structural  architecture  of  the  breast.  The  spatial  distribution  and  the  relative  concentrations 
of  these  chromophores  are  patient  dependent  and  even  more  can  vary  over  time  in  the  same 
patient  due  to  hormonal  regulation  (Chance  2001,  Cubeddu  et  al  2000,  Durduran  et  al  2002, 
Shah  et  al  2001,  Srinivasan  et  al  2003). 

To  enhance  diffuse  optical  tomography  (DOT)  performances,  researchers  have  proposed 
fusing  optical  techniques  with  other  medical  imaging  modality.  Magnetic  resonance  imaging 
(MRI)  is  the  perfect  candidate  for  optical  co-registration  (Brooksby  etal  2003,  Guven  etal 
2004,  Ntziachristos  et  al  2002,  Pei  et  al  1999,  Pogue  and  Paulsen  1998).  MRI  provides  high 
spatial  resolution  maps  of  the  breast  optical  structure  that  are  relevant  to  the  water  and  lipid 
distribution.  Moreover,  MRI  can  provide  a  means  to  estimate  the  concentration  of  these  two 
structural  chromophores  (Merrit  etal  2003).  In  this  paper,  we  investigate  the  first  step  towards 
incorporating  physiological  and  spatial  a  priori  information  derived  from  MRI. 


2.  Methods 


2.1 .  Forward  model 


The  propagation  of  NIR  light  in  tissue  is  well  modelled  by  the  diffusion  equation.  In  the  case 
of  heterogeneity,  the  diffusion  equation  can  be  solved  by  a  perturbative  approach  (O’Leary 
1996).  In  this  work,  we  have  used  the  Rytov  approximation  approach.  In  the  case  of  DOT, 
multiple  source-detector  pairs  are  used.  The  medium  under  consideration  is  sampled  in  voxels 
and  the  problem  can  be  written  as  a  matrix  equation  (Arridge  1999),  i.e.: 
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where  $  ( rsi ,  rdi )  is  the  diffuse  perturbative  phase  for  the  ith  source-detector  pair,  Wjj  (O  ’Leary 
1996)  is  the  weight  for  the  yth  voxel  and  the  ith  source-detector  pair  and  <S/ra(r;  )  is  the 
differential  absorption  coefficient  of  the  jth  voxel.  We  limited  our  problem  to  image  the 
absorption  coefficient.  Boundary  conditions  for  semi-infinite  geometries  and  slab  geometries 
are  derived  using  the  extrapolated  boundary  condition  and  the  image  source  technique 
(Haskell  et  al  1994). 


2.2.  Functional  imaging 

The  estimation  of  the  absorption  coefficient  at  several  wavelengths  enables  the  provision  of 
spatial  maps  of  the  targeted  chromophores.  In  the  case  of  the  breast,  the  four  chromophores 
of  potential  diagnostic  interest  are  the  oxy-  and  deoxy-haemoglobin,  the  water  and  the  lipid. 
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The  concentrations  of  these  breast  constituents  are  linearly  related  to  the  absorption  values 
through  the  linear  system: 
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where  8ixkk(rj)  is  the  differential  absorption  coefficient  at  the  &th  wavelength  Xk  and  for  the 
;th  voxel,  s^k  is  the  extinction  coefficient  of  the  Cth  chromophore  at  the  Ath  wavelength  and 
<5[C](rj)  is  the  differential  concentration  of  the  Cth  chromophore  for  the  jth  voxel. 

Solving  this  linear  system  on  a  voxel  basis  for  a  spectral  set  of  absorption  distributions 
gives  the  functional  maps  required  for  diagnostic  purposes.  We  will  refer  to  this  approach  as 
‘indirect  imaging’  through  this  paper. 

Recently,  to  overcome  the  burden  and  minimize  the  systematic  errors  due  to  the 
ill-condition  of  both  inverse  problems,  a  method  aiming  at  directly  imaging  the  functional 
parameters  has  been  proposed  (Durduran  et  al  2001).  This  method  takes  advantage  of  the 
linear  relationship  described  previously  and  formulates  the  inverse  problem  as 
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This  new  linear  system  is  poorly  conditioned  and  great  care  should  be  taken  during  the 
pre-conditioning  of  this  sensitivity  matrix.  In  our  case,  we  used  an  average  column  scheme 
(Pei  etal  2001)  applied  to  each  sub-matrix  of  the  kernel  as  written  in  (3).  This  specific  pre¬ 
conditioning  scheme  led  to  the  most  accurate  reconstructions.  We  will  refer  to  the  scheme 
expressed  in  (3)  as  ‘direct  imaging’  through  this  document. 

The  functional  parameters  that  we  are  interested  in  are  the  blood  volume: 

[BV]  =  [Hb]  +  [Hb02]  (4) 

and  the  relative  saturation: 

[Hb02] 


[Sa02]  = 


[Hb]  +  [Hb02] ' 


(5) 


2.3.  Bayesian  framework 

Due  to  the  ill-posed  and/or  underdetermined  nature  of  the  DOT  problem,  the  solution  of  the 
inverse  problem  is  not  typically  robust.  One  avenue  to  overcome  this  difficulty  is  to  incorporate 
a  priori  information  constraining  the  space  of  unknowns. 

The  Bayesian  approach  provides  a  natural  framework  to  incorporate  prior  information. 
Guven  et  al  (2004)  proposed  an  algorithm  based  on  the  Bayesian  framework  with  a  spatially 
varying  a  priori  probability  density  function  extracted  from  MRI  anatomical  maps.  Here  we 
propose  to  extend  this  algorithm  to  a  spatial  physiological  prior.  The  full  derivation  of  the 
theoretical  developments  for  the  algorithm  is  presented  in  detail  in  Guven  et  al  (2004)  and  we 
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follow  herein  the  same  mathematical  expressions.  We  present  in  this  investigation  only  the 
salient  features  of  the  theoretical  approach  that  are  relevant  to  this  work. 

The  available  high-resolution  anatomical  image  is  segmented  into  sub-images  that 
represent  major  tissue  types  (typically:  parenchyma,  glandular  and  tumour).  Prior  probability 
density  function  of  the  image  is  formulated  in  such  a  way  that  each  sub-image  is  assigned  a 
mean  value  that  need  not  be  equal  to  its  actual  optical  value;  and  a  ‘confidence  level’  is  defined 
in  the  form  of  an  image  variance  formulation  to  allow  local  variations  within  sub-images.  As 
a  consequence,  the  overall  formulation  of  the  prior  information  becomes  spatially  varying, 
which  is  specific  to  the  image  of  interest.  Maximum  a  posteriori  (MAP)  estimate  of  the  image 
is  formed  based  on  the  formulation  of  the  image’s  probability  density  function 

xMap  =  arg  maxflog  p(y  |x)  +  log  p(x)}  (6) 

X 

where  p(y\x)  is  the  data  likelihood  function  and  p(x)  is  the  probability  density  function  of  the 
corresponding  image.  An  ‘alternating  minimization’  algorithm,  which  sequentially  updates 
the  unknown  parameters,  is  used  to  solve  the  resulting  optimization  problem. 

For  our  purposes,  the  probability  density  function  of  the  fth  sub-image,  as  defined  in  the 
spatial  prior,  becomes 

p(xi  |ff,)  =  - - \-jTii  exP  (~  Hx«'  "  ci  II2)  i  =  1,2,...,  M  (7) 

(2jrcr?)  l/  \  2of  ) 

where  M  is  the  number  of  sub-regions  and  Ni  is  the  number  of  voxels  in  the  ith  sub-image, 
Xi  is  the  unknown  sub-image,  C,-  is  the  assigned  chromophore  mean  concentration  and  of 
the  single  variance,  which  is  an  unknown  parameter  (estimated  during  the  solution).  To 
incorporate  the  confidence  level  into  the  statistical  reconstruction  procedure,  the  sub-image 
variances  are  expressed  as 

Pfo)  =  {2nyf)m  CXP  (~ ll<r‘  ""  &i l|2)  i  —  \,2, . . . ,  M  (8) 

where  y,  is  the  variance  and  ff,-  the  mean  value  of  ff,-.  These  two  last  parameters  are  a  priori 
defined  by  the  user.  Thus  the  confidence  level  incorporated  into  the  statistical  reconstruction 
procedure  is  defined  on  the  chromophores  concentration.  Hence,  physiological  priors  are 
implicitly  defined  and  allow  constraining  the  reconstruction  in  physiologically  meaningful 
ranges. 

2.4.  Measurement  generation 

Measurements  were  obtained  by  solving  the  frequency-domain  diffusion  equation  with  a  finite 
difference  approach  (FDM).  We  restricted  our  simulations  to  a  two-dimensional  (2D)  geometry 
for  computational  efficiency  and  to  a  continuous  wave  data  set  type.  The  slab  thickness  was 
6  cm  simulating  a  soft  compressed  breast.  We  placed  nine  sources  on  one  side  of  the  slab  and 
nine  detectors  on  the  other  side,  both  evenly  stretched  along  8  cm  (cf  figure  1).  Two  square 
inclusions  of  1  cm2  were  simulated  in  this  model. 

The  functional  properties  were  chosen  to  mimic  typical  values  encountered  in  the  human 
breast.  Table  1  summarizes  the  functional  parameters  simulated. 

The  optical  models  were  computed  for  six  wavelengths  to  replicate  the  spectral 
information  gathered  by  our  time  resolved  instrument  that  performs  co-registration  with  MRI 
(Intes  et  al  2002).  The  subsequent  optical  properties  are  compiled  in  table  2.  The  FDM 
computations  were  performed  with  a  1  mm  mesh  resolution.  The  sources  and  detectors  were 
positioned  2  cm  away  from  the  edges  to  avoid  any  boundary  effects. 
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Figure  1.  Optical  model  used  for  the  simulations.  The  optical  properties  displayed  here  correspond 
to  X  —  750  nm. 


Table  1.  Functional  parameters  used  to  define  the  optical  properties  to  generate  the  synthetic 
measurements. 


Background 

Left  object 

Right  object 

BV  (mM) 

0.02 

0.06 

0.12 

Sa02  (%) 

70 

65 

50 

Table  2.  Optical  parameters  for  the  spectral  set  investigated  herein.  No  [H2O]  and  [Li]  constituents 
were  used  in  these  simulations  for  simplicity. 


690 

750 

780 

805 

830 

850 

Ms  (cm-1) 

9.268 

8.386 

8.000 

7.646 

7.425 

7.216 

Makg  (cm-1) 

0.016 

0.016 

0.016 

0.016 

0.018 

0.019 

M^ft  (cm"1) 

0.054 

0.050 

0.050 

0.049 

0.053 

0.056 

p,f'8h‘  (cm-1) 

0.140 

0.115 

0.107 

0.095 

0.100 

0.105 

3.  Results 

First,  the  estimation  of  the  functional  parameters  simulated  is  accomplished  with  a  conjugate 
gradient  descent  (CGD)  algorithm  in  the  case  of  indirect  imaging  without  using  any  kind  of 
prior.  The  results  are  displayed  in  figure  2. 

Then  reconstructions  using  the  indirect  imaging  approach  within  the  Bayesian  framework 
are  provided  in  figure  3 .  First,  an  estimation  of  the  absorption  coefficients  at  the  six  wavelengths 
was  performed  using  spatial  a  priori  information  and  a  conjugate  gradient  algorithm  with  the 
Polak-Ribiere  method  (Polak  and  Ribiere  1969).  The  absorption  coefficients  mean  values 
assigned  were  the  simulated  ones  with  a  30%  level  of  confidence  (Guven  et  al  2004).  Then 
classical  spectroscopy  was  performed  on  the  resulting  optical  maps  using  the  linear  system 
of  (3)  reduced  to  [Hb02]  and  [Hb]  concentrations.  Last,  the  results  using  physiological  and 
spatial  a  priori  information  are  depicted  in  figure  4.  The  chromophore  concentration  means 
assigned  were  the  exact  ones  with  a  30%  level  of  confidence. 

The  functional  quantitative  values  retrieved  from  each  case  and  for  the  relevant  region  of 
interest  are  summarized  in  table  3. 
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Figure  2.  (a)  Blood  volume  (mM)  reconstructions  and  (b)  saturation  estimates  (%)  in  the  case  of 
classical  indirect  imaging. 
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Figure  3.  (a)  Blood  volume  (mM)  reconstructions  and  (b)  saturation  estimates  (%)  in  the  case 
of  Bayesian  indirect  imaging.  The  assigned  mean  values  of  the  absorption  correspond  to  the  true 
absorption  values  with  a  30%  level  of  confidence. 


Figure  4.  (a)  Blood  volume  (mM)  reconstructions  and  (b)  saturation  estimates  (%)  in  the  case  of 
Bayesian  direct  imaging.  The  assigned  mean  values  correspond  to  the  true  concentration  values 
with  a  30%  level  of  confidence. 


4.  Discussion 

In  all  the  cases  investigated  herein,  the  two  inclusions  were  reconstructed.  As  expected, 
the  incorporation  of  spatial  a  priori  information  increased  the  image  resolution.  In  the 
case  of  indirect  imaging  without  a  priori  information  the  two  inclusions  are  elongated  and 
thus  the  contrast  is  diluted.  This  fact  is  reflected  in  the  estimation  of  the  blood  volume  that  is 
underestimated  for  both  objects  (~35%  lower  estimate).  When  a  priori  anatomical  information 
is  incorporated  within  the  Bayesian  framework,  the  estimation  of  the  blood  volume  is  more 
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Table  3.  Functional  parameters  recovered  with  the  three  approaches  and  for  the  three  different 
functional  areas.  The  values  proposed  here  correspond  to  the  mean  value  of  the  entire  ROI  defined 
as  the  a  priori  spatial  masks.  [BV]  in  mM  and  [Sa02]  in  %. 


Background 

Left  object 

Right  object 

[BV] 

[Sa02] 

[BV] 

[Sa02] 

[BV] 

[Sa02] 

Indirect  imaging 

0.023 

71 

0.042 

67 

0.075 

57 

Bayesian  indirect  imaging 

0.019 

72.46 

0.065 

63.73 

0.122 

55.44 

Bayesian  direct  imaging 

0.020 

69.95 

0.061 

65.27 

0.127 

50.02 

accurate  falling  to  less  than  10%  of  misestimation  in  the  worse  case  (1%  in  the  best  case). 
Moreover,  when  direct  imaging  is  applied  in  this  framework,  the  two  objects  are  recovered 
respectively  within  1%  and  5%  for  the  left  and  right  objects  compared  to  8%  and  2%  in  the 
case  of  indirect  imaging. 

The  enhancement  of  the  technique  is  even  more  notable  when  SaC>2  is  considered.  The 
recovery  of  SaC>2  is  then  more  robust  but  also  less  prone  to  artefacts.  Especially  in  the  case  of 
the  right  object,  the  recovery  of  the  SaC>2  is  superior  to  the  direct  imaging  approach  within  the 
Bayesian  framework. 

Overall,  these  simulations  highlight  the  potential  of  our  approach  to  provide  more  accurate 
maps  of  the  relevant  functional  diagnostic  parameters.  These  simulations  were  limited  to 
[BV]  and  Sa02  recovery.  The  approach  is  easily  extensible  to  [H20]  and  [Li].  Especially,  the 
ability  of  MRI  to  provide  spatial  concentration  of  these  chromophores  for  the  assigned  prior 
strengthens  our  approach  (Merrit  et  al  2003). 

The  mean  values  of  the  chromophores  and  the  anatomical  a  priori  information  are  defined 
by  the  user.  The  anatomical  maps  are  provided  by  the  structural  MRI  maps  (such  as  Tl) 
and  tumour  delineation  could  be  performed  by  Gd  enhanced  MRI.  Moreover,  estimation  of 
the  water,  lipid  and  BV  concentrations  are  feasible  with  MRI.  These  values  can  be  used  for 
the  mean  concentration  priors.  Also,  an  alternative  is  to  obtain  these  values  with  a  simple 
algorithm  such  as  diffuse  optical  spectroscopy  (Ntziachristos  et  al  2002).  The  paradigm  of 
uniform  concentrations  for  a  certain  tissue  type  is  overcome  in  the  Bayesian  formulation  by 
defining  a  level  of  confidence.  These  effects  of  confidence  ensure  the  reconstructions  are  in  a 
physiologically  meaningful  range  and  allow  the  recovery  of  heterogeneous  concentration  in  a 
certain  kind  of  tissue  type.  Such  heterogeneous  estimates  are  visible  in  figures  3  and  4. 

This  last  point  is  of  paramount  importance  in  the  case  of  MRI-assisted  DOT.  The  contrast 
provided  by  MRI  is  not  expected  to  be  exactly  congruent  with  the  optical  contrast.  Allowing 
heterogeneous  bounded  reconstructions  within  each  sub-image  compensates  for  discrepancies 
between  both  modalities. 

5.  Conclusion 

We  reported  in  this  work  our  first  step  towards  incorporating  physiological  and  spatial 
a  priori  information  derived  from  MRI  to  assist  DOT.  Better  estimates  of  the  main  functional 
parameters  that  are  [BV]  and  [Sa02]  were  achieved.  More  accurate  functional  maps  can  lead 
to  an  increase  of  the  sensitivity  and  specificity  of  optical  techniques  by  incorporating  structural 
contrasts  that  are  known  to  occur  in  breast  cancer.  Moreover,  these  accurate  estimations  will 
be  important  for  application  such  as  therapy  monitoring  (Zhang  et  al  2003). 

This  preliminary  work  will  be  continued  with  experimental  validation  and  incorporated 
in  our  ongoing  clinical  trial  at  the  University  of  Pennsylvania.  Also,  the  investigation  will  be 
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extended  to  the  potential  of  our  approach  to  recover  exogenous  contrast  agent  concentrations. 
New  classes  of  contrast  agents  that  are  tumour  specific  (Licha  2002)  and/or  compatible  with 
both  technologies  (Josephon  et  al  2002)  will  be  investigated.  The  inclusion  of  spatial  and 
temporal  a  priori  information  is  expected  to  provide  a  better  estimate  of  the  fluorochrome 
concentration  but  also  to  probe  physiological  parameters  such  as  extravasation  (Intes  et  al 
2003,  Cuccia  et  al  2003)  that  are  still  elusive  with  DOT  for  thick  tissue  investigation. 
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Abstract 


Compartmental  modeling  of  indocyanine  green  (ICG)  pharmacokinetics,  as  measured  by  near 
infrared  (NIR)  techniques,  has  the  potential  to  provide  diagnostic  information  for  tumor  dif¬ 
ferentiation.  In  this  paper,  we  present  three  different  compartmental  models  to  model  the 

pharmacokinetics  of  ICG  in  cancerous  tumors.  We  introduce  a  systematic  and  robust  ap- 

i 

proach  to  model  and  analyze  ICG  pharmacokinetics  based  on  the  extended  Kalman  filtering 

| 

(EKF)  framework.  The  proposed  EKF  framework  effectively  models  multiple-compartment 
and  multiple-measurement  systems  in  the  presence  of  measurement  noise  and  uncertainties 
in  model  dynamics.  It  provides  simultaneous  estimation  of  pharmacokinetic  parameters  and 
ICG  concentrations  in  each  compartment.  Moreover,  the  recursive  nature  of  the  Kalman  filter 

estimator  potentially  allows  real  time  monitoring  of  time  varying  pharmacokinetic  rates  and 

i 

concentration  changes  in  different  compartments.  Additionally,  we  introduce  an  information 
theoretic  criteria  for  the  best  compartmental  model  order  selection,  and  residual  analysis  for 
the  statistical  validation  of  the  estimates.  We  tested  our  approach  using  the  ICG  concentration 
data  acquired  from  four  Fischer  rats  carrying  adenocarcinoma  tumor  cells.  Our  study  indicates 

I 

I 

that,  in  addition  to  the  pharmacokinetic  rates,  EKF  model  may  provide  parameters  that  may  be 
useful  for  tumor  differentiation. 

Index  Terms 

Extended  Kalman  filter,  indocyanine  green,  compartmental  analysis,  pharmacokinetics,  tumor 
characterization. 
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I.  Introduction 


Near  infrared  (NIR)  diffuse  optical  imaging  and  spectroscopy  methods  provide  quantitative  functional 
information  that  cannot  be  obtained  by  the  conventional  radiological  methods  [1],  [2],  [3].  NIR  techniques 
can  provide  in  vivo  measurements  of  the  oxygenation  and  vascularization  states,  uptake  and  release 
of  optical  contrast  agents,  and  chromophore  concentrations  with  high  sensitivity.  In  particular,  NIR 
diffuse  optical  techniques  in  conjunction  with  optical  contrast  agents  have  the  potential  to  characterize 
angiogenesis,  and  to  differentiate  between  malignant  and  benign  tumors  [4],  [5],  [6],  [7]. 

At  present,  indocyanine  green  (ICG)  is  the  only  NIR  optical  agent  approved  for  human  use.  In  NIR 
measurements,  the  presence  of  ICG  within  an  imaging  volume  results  in  an  increased  signal  that  can  be 
observed  over  the  course  of  the  experiment.  Study  of  the  time  kinetics  of  ICG  concentration  curves  may 
provide  physiologically  relevant  information  for  tumor  differentiation.  Specifically,  cancerous  tissue  types 
are  expected  to  show  high  and  fast  uptake  due  to  the  proliferation  of  ’’leaky”  angiogenetic  microvessels, 
while  normal  and  fatty  tissue  show  little  uptake. 

A  number  of  research  groups  reported  compartmental  modeling  of  ICG  time-kinetic  measurements 
using  NIR  methods  for  tumor  diagnosis  in  animal  and  human  subjects  [8],  [9],  [10].  Compartmental 
model  is  a  mathematical  description  of  the  concentrations  of  contrast  agents  in  which  each  compartment 
represents  kinetically  distinct  tissue  type.  It  consists  of  a  set  of  coupled  ordinary  partial  differential 
equations  (ODE)  and  a  measurement  model.  Coefficients  of  the  ODE’s  are  the  physiological  parameters 
of  interest  that  represent  rates  of  exchange  between  different  compartments.  These  parameters  are  non- 
linearly  related  to  the  total  concentration  of  ICG  measured  by  NIR  methods.  Furthermore,  concentration 
of  ICG  in  each  compartment  cannot  be  directly  measured  non  invasively  by  NIR  techniques,  making  the 
pharmacokinetic  parameter  estimation  a  highly  non-linear  problem. 

Current  methods  of  ICG  compartmental  modeling  involve  curve  fitting  methods  and  various  techniques 
for  solving  differential  equations.  Gurfinkel  et  al.  presented  a  two-compartment  model  for  ICG  kinetics 
and  estimated  model  parameters  [8].  The  measurements  were  obtained  using  a  frequency  domain  photon 
migration  system  coupled  with  a  charge-coupled  device.  The  pharmacokinetic  parameters  are  estimated 
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for  each  pixel  based  on  a  curve  fitting  method.  This  study  indicates  that  model  parameters  show  no 
difference  in  the  ICG  uptake  rates  between  normal  and  diseased  tissue.  Cuccia  et  al.  presented  a  study  of 
the  dynamics  of  ICG  in  an  adenocarcinoma  rat  tumor  model  [9].  A  two-compartmental  model  describing 
the  ICG  dynamics  is  used  to  quantify  physiologic  parameters  related  to  capillary  permeability.  The  ICG 
concentration  curves  were  fitted  to  the  compartmental  model  using  a  non-linear  least  squares  Levenberg- 
Marquart  algorithm.  It  was  shown  that  different  tumor  types  have  different  capillary  permeability  rates. 
Intes  et  al.  presented  the  uptake  of  ICG  by  breast  tumors  using  a  continuous  wave  diffuse  optical 
tomography  apparatus  [10].  A  two-compartment  model  is  used  to  analyze  the  pharmacokinetics  of  ICG. 
A  curve  fitting  algorithm,  namely  the  non-linear  Nelder-Mead  simplex  search,  is  used  to  estimate  the 
pharmacokinetic  parameters.  This  study  shows  that  the  malignant  cases  exhibit  slower  rate  constants 
(uptake  and  outflow)  as  compared  to  healthy  tissue. 

While  the  studies  described  above  demonstrate  the  feasibility  of  the  ICG  pharmacokinetics  in  tumor 
characterization;  due  to  highly  non-linear  nature  of  the  pharmacokinetic  parameter  estimation,  variation  in 
parameter  values  from  one  subject  to  another,  and  sparse  data  available  in  clinical  and  laboratory  settings, 
a  systematic  and  robust  approach  is  needed  to  model,  estimate  and  analyze  ICG  pharmacokinetics.  Such 
an  approach  must  include:  i )  a  method  for  compartmental  model  order  selection,  ii)  a  robust  method  of 
estimating  ICG  pharmacokinetic  parameters,  and  Hi)  a  method  of  validating  the  selected  model  and  the 
estimation  results. 

In  this  paper,  we  present  three  different  compartmental  models  for  the  ICG  pharmacokinetics  in 
cancerous  tumors  and  propose  an  extended  Kalman  filtering  (EKF)  framework  to  estimate  the  model 
parameters.  The  models  capture  the  transportation  of  ICG  between  the  vascular  and  extravascular  com¬ 
partments,  including  interstitial  fluid  region,  parenchymal  cell,  intracellular  binding  site,  and  extravascular, 
extracellular  spaces  (EES).  Kalman  filter  (KF)  is  an  optimal  recursive  modeling  and  estimation  method 
with  numerous  advantages  in  ICG  pharmacokinetic  modeling.  These  include:  i )  effective  modeling  of 
multiple  compartments,  and  multiple  measurement  systems  governed  by  coupled  ordinary  differential 
equations,  in  the  presence  of  measurement  noise  and  uncertainties  in  the  compartmental  model  dynam- 
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ics;  ii )  simultaneous  estimation  of  pharmacokinetic  model  parameters  and  ICG  concentrations  in  each 
compartment,  which  is  not  accessible  in  vivo  by  means  of  NIR  techniques;  iii)  recursive  estimation  of 
time-varying  pharmacokinetic  model  parameters;  iv)  statistical  validation  of  estimated  concentrations  and 
error  bounds  on  the  pharmacokinetic  parameter  estimates;  v)  incorporation  of  available  a  priori  information 
about  the  initial  conditions  of  the  permeability  rates  into  the  estimation  procedure;  vi )  potential  real-time 
monitoring  of  ICG  pharmacokinetic  parameters  and  ICG  concentrations  in  different  compartments  due 
to  the  recursive  nature  of  EKF  estimation  method.  Additionally,  we  present  a  method  of  selecting  the 
optimal  compartmental  model  order  based  on  Bayesian  information  criteria,  and  a  statistical  validation 
method  based  on  the  residual  analysis. 

We  tested  our  approach  using  the  ICG  concentration  data  acquired  from  four  Fisher  rats  carrying  adeno¬ 
carcinoma  tumor  cells.  Two-,  three-  and  four-compartment  models  are  fitted  to  data  and  pharmacokinetic 
model  parameters  and  concentrations  in  different  compartments  are  estimated  using  the  EKF  framework. 
Bayesian  information  criteria  suggests  that  the  two-compartment  model  provides  a  sufficient  fit  for  our 
data.  The  estimated  model  order  and  the  model  parameters  are  further  validated  by  the  residual  analysis. 
The  model  parameters  were  used  to  differentiate  between  two  types  of  cancerous  tumors.  Our  study 
suggests  that  the  permeability  rates  out  of  the  vasculature  are  higher  in  edematous  tumors  as  compared 
to  necrotic  tumors.  Additionally,  we  observed  that  in  the  two-compartment  model,  the  ICG  concentration 
curve  is  higher  in  the  EES  compartment  in  edematous  tumors.  This  suggests  that  the  ratio  of  the  peak 
value  of  the  ICG  concentrations  in  different  compartments  may  be  a  useful  parameter  to  differentiate 
tumors. 

The  paper  is  organized  as  follows:  In  Section  II,  we  present  the  two-,  three-  and  four-compartment 
models  for  the  ICG  pharmacokinetics  in  tissue.  In  Section  III,  we  present  the  state-space  representation  of 
the  compartmental  models;  estimation  of  ICG  pharmacokinetics  parameters  and  ICG  concentrations  in  the 
EKF  framework;  and  optimal  model  order  selection  criteria.  In  Section  IV,  we  present  the  experimental 
results  obtained  from  Fischer  rat  data.  Section  V  summarizes  our  results  and  conclusion.  The  appendix 
includes  the  derivation  of  the  likelihood  function  used  in  the  Bayesian  information  criteria. 
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II.  ICG  Pharmacokinetic  Modeling  Using  NIR  MEASUREMENTS 

A.  Indocyanine  Green 

ICG  is  an  optical  dye  commonly  used  in  retinopathy  and  hepatic  diagnostics.  Given  its  low  toxicity 
and  FDA  approval,  it  has  recently  been  utilized  as  a  blood  pooling  agent  for  the  detection  and  diagnosis 
of  cancerous  tumors  by  means  of  NIR  optical  methods.  The  absorption  peak  of  ICG  is  805  nm  and 
the  fluorescence  peak  is  at  830  nm.  ICG  has  strong  affinity  for  blood  proteins.  In  plasma,  ICG  is  near- 
completely  bound,  primarily  to  albumin.  As  a  result,  its  in  vivo  kinetics  are  similar  to  those  of  a  70  kD 
molecule,  although  it  has  a  molecular  weight  of  about  700  D  [11],  [12],  [13],  [14],  [15]. 

ICG  is  eliminated  from  the  body  primarily  through  the  bile.  Outside  of  the  circulatory  system,  it  is  not 
available  for  removal  until  it  returns  to  the  system.  The  kinetics  of  this  transition  offers  a  potential  means 
of  non-invasively  assessing  the  leakiness  of  large  molecules  from  the  microvasculature;  this  permeability 
is  a  characteristic  of  the  poorly  developed  vasculature  observed  in  angiogenesis.  The  increase  in  local 
microvasculature  density  is  also  expected  to  induce  increased  perturbation  in  the  optical  signal  from 
intercapillary  ICG. 

There  are  some  differences  in  the  delivery  of  ICG  between  normal  and  cancerous  vasculature.  In  normal 
tissue,  ICG  acts  as  a  blood  flow  indicator  in  tight  capillaries  of  normal  vessels.  However  in  tumor,  ICG 
may  act  as  a  diffusible  (extravascular)  flow  in  the  leaky  capillary  of  cancer  vessels.  Additionally,  the 
permeability  rate  is  expected  to  increase  as  the  malignancy  advances  [9],  [10],  Figure  1  (a)  and  (b) 
illustrates  the  ICG  flow  for  healthy  and  malignant  tissue,  respectively. 

B.  Compartmental  Analysis  of  ICG  Pharmacokinetics 

Compartmental  modeling  allows  relatively  simple  and  effective  mathematical  representation  of  complex 
biological  responses  due  to  contrast  agents.  A  region  of  interest  is  assumed  to  consist  of  a  number  of 
compartments,  generally  representing  a  volume  or  a  group  of  similar  tissues  into  which  the  contrast 
agent  is  distributed.  The  concentration  change  in  a  specific  compartment  is  modeled  as  a  result  of  the 
exchange  of  contrast  agent  between  connected  compartments.  These  changes  are  modeled  by  a  collection 
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of  coupled  ODEs;  each  equation  describing  the  time  change  dictated  by  the  biological  laws  that  govern 
the  concentration  exchanges  between  the  interacting  compartments  [16],  [17],  [18],  [19].  In  this  work,  we 
investigate  three  different  compartmental  models  for  the  ICG  kinetics  and  determine  the  optimal  model 
order,  based  on  Bayesian  information  criteria. 

1)  The  four-compartment  model:  Figure  2  illustrates  the  capillary  and  extracapillary  space  relevant 
to  the  four  compartment  model.  The  four-compartment  model  includes  capillary  region,  interstitial  fluid 
region,  parenchymal  cell  region  and  intracellular  binding  site  as  compartments  [20],  The  ICG,  injected 
intravenously  into  the  subject,  can  pass  through  from  capillary  into  reversible  binding  site  inside  the 
cell  through  the  interstitial  fluid  region  and  the  parenchymal  cell  region  [20],  [21],  [22],  Moreover, 
in  advanced  tumor  stages,  the  leakiness  around  the  tumor  vessels  is  expected  to  increase,  resulting  in 
higher  permeability  rates  during  the  transportation  of  ICG  into  the  compartments.  A  block  diagram  of 
the  four-compartment  transport  and  chemical  model  of  ICG  delivery  is  shown  in  Figure  3(a). 

Let  Cp,  Ci,  Cpc,  Cb  denote  the  ICG  concentrations  in  plasma,  interstitial  fluid  region,  parenchymal 
cell  region  and  intracellular  binding  site,  respectively;  and  let  k$t,  ^  ^  ,  k^\  ki^  and 

be  the  constants  used  as  equilibrium  coefficients  as  shown  in  Figure  3(a).  Then  the  set  of  differential 
equations  representing  the  ICG  transition  between  the  four  compartments  are  given  as  follows: 

The  leakage  into  and  the  drainage  out  of  plasma: 

=  k^Cft)  -  k^Cp(t)  -  k(04JtCp(t).  (1) 

The  leakage  into  and  the  drainage  out  of  interstitial  fluid  region: 

—fp-  =  k{4)Cp(t)  -  k{4)Ci{t )  -  k^Cft)  +  kfcpc{t).  (2) 

The  leakage  into  and  the  drainage  out  of  parenchymal  cell: 

=  k^Ci(t)  -  kfcpc{t)  -  k^Cpcit)  +  kfcb{t).  (3) 

The  leakage  into  and  the  drainage  out  of  intracellular  binding  site: 

=  k^Cpc(t)  -  kfcb{t).  (4) 
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Physiologically,  the  equilibrium  constants  are  defined  by  the  permeability  surface  area  products  given  as 
PSp,  where  P  is  the  capillary  permeability  constant,  S  is  the  capillary  surface  area,  and  p  is  the  tissue 
density.  k^Jt  is  proportional  to  the  flow  rate  into  and  out  of  the  capillary  and  k^\  kb4\  ktf*  ,  k(f\  ki4\ 
and  k^p  represent  intra-tissue  physiologic  effects  during  ICG  delivery  from  capillary  to  binding  site.  Note 
that  the  superscript  denotes  the  order  of  the  compartmental  model. 

The  actual  bulk  ICG  concentration  in  the  tissue  measured  by  NIR  spectroscopy,  m(i),  is  a  linear 
combination  of  the  ICG  concentrations  in  four  different  compartments. 

m(t)  =  v^Cp(t)  +  v^Ci(t)  +  vj$Cpe(t)  +  v(b4)Cb(t),  (5) 

where  Vp4\  v-4\  vp4\  vb4\  are  volume  fractions  of  plasma,  interstitial  fluid  region,  parenchymal  cell 
region  and  intracellular  binding  site,  respectively. 

2)  The  three-compartment  model:  In  this  model,  the  parenchymal  cell  and  intracellular  binding  site 
compartments  are  combined  to  form  a  single  compartment  called  parenchymal  cell.  This  amounts  to  the 
assumptions  that  the  transition  of  ICG  into  the  intracellular  binding  site  is  negligible  as  compared  to  the 
other  compartments,  and  therefore  omitted  from  the  model.  A  block  diagram  of  the  three-compartment 
transport  and  chemical  model  of  ICG  delivery  is  shown  in  Figure  3(b).  Three-compartment  transport 
equations  are  given  as  follows: 

The  leakage  into  and  the  drainage  out  of  plasma: 

=  kf] Ci.it)  -  k^Cpit)  -  k^Cpit)  (6) 

The  leakage  into  and  the  drainage  out  of  interstitial  space: 

-  k^Cpit)  -  fcfa(f)  -  k^Ci(t)  +  kfcpc{t )  (7) 

The  leakage  into  and  the  drainage  out  of  parenchymal  cell: 

=  kjPCiit)  -  kfcpc{t)  (8) 

The  total  ICG  concentration  measured  by  NIR: 

m(t)  =  v^Cp(t)  +  vf]Ci{t )  +  v$Cpc(t)  (9) 

where  \  v(f  \  and  Cp,  C\,  Cpc  are  as  defined  in  the  four-compartment  model. 
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3)  The  two-compartment  model:  In  the  two-compartment  model,  the  tumor  region  is  assumed  to  be 
composed  of  two  compartments,  namely  the  plasma  and  EES  [9],  [23],  [24].  EES  is  defined  as  the  region 
that  lies  outside  of  both  the  vascular  region  and  the  tumor  cells.  The  transition  of  the  ICG  to  the  third 
and  fourth  compartments  are  assumed  to  be  negligible.  Therefore  the  last  two  compartments  in  the  four 
compartment  model  is  omitted.  We  consider  transcapillary  leakage  to  occur  only  at  the  tumor  site.  We 
also  assume  that  a  small  perturbation  of  the  global  plasma  concentration  does  not  affect  the  bulk  removal. 
Figure  3(c)  shows  the  block  diagram  of  the  two-compartment  model  for  the  ICG  kinetics.  Let  Cp  and 
Ce  denote  the  ICG  concentrations  in  the  plasma  and  EES,  respectively.  Then  the  two-compartment  ICG 
chemical  transport  equations  are  given  as  follows: 

The  leakage  into  and  the  drainage  out  of  plasma: 

<^fL  =  k^Ce(t)  -  k^Cp(t)  -  k™cp(t).  (10) 

The  leakage  into  and  the  drainage  out  of  the  EES: 

^M  =  kPcp(t)-k^ce(t).  (ii) 

The  parameters  k^  and  k^  govern  the  leakage  into  and  the  drainage  out  of  the  EES,  respectively. 
The  parameter  k^Jt  describes  the  ICG  elimination  from  the  body  through  kidneys  and  liver. 

Actual  bulk  ICG  concentration  in  the  tissue  measured  by  NIR  is  a  linear  combination  of  plasma  and 
EES  ICG  concentrations  given  by: 

m(t)  =  42)Cp(f)  +  v^Ceit),  (12) 

where  the  parameters  and  vj>2)  denote  the  plasma  and  EES  volume  fractions,  respectively. 

III.  Extended  Kalman  Filtering  for  the  ICG  Pharmacokinetics 

For  the  rest  of  our  discussion,  we  shall  use  the  explicit  form  of  the  two-compartment  model  as  a 
running  example  to  clarify  our  notation. 
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A.  State-space  Representation  of  the  ICG  Pharmacokinetics 

Coupled  differential  equations  resulting  from  the  two-compartment  modeling  of  the  ICG  pharmacoki¬ 
netics  can  be  expressed  in  state-space  representation  as  follows: 


Ce(t) 

u(2) 

~Kb 

- - 

S" 

se. 

Ce(t) 

_  cP(t)  _ 

jl(2) 

Kb 

-(kP  +  fcS)  _ 

_  cp(t)  _ 

(13) 


m(t) 


42)  42) 


Ce(t) 

CP(t) 


+  v(t) 


where  u>(t)  and  p(t)  are  uncorrelated  zero  mean  Gaussian  processes  with  covariance  matrix  Q,  and 
variance  cr2,  respectively. 

The  closed  form  of  the  continuous  time  state-space  representation  for  the  n— compartment  model  is 
given  by: 

dC(t )  =  K(an)C(t)dt  +  u>(t)dt, 


m(t )  =  V(an)C(t)  +  r/(t). 


(14) 


In  equation  (14),  C(f)  denotes  the  concentration  vector;  K(an)  is  the  KF  system  matrix,  V(an)  is 
the  KF  measurement  matrix  and  an  is  the  parameter  vector  whose  elements  are  the  pharmacokinetic 
constants  and  the  volume  fractions  for  the  n— compartment  model.  For  example  the  parameter  vector  a.2 
for  the  two-compartment  model  is  given  as 

«2  =  [*42)  fcfe2)  V(e}  42)].  (15) 

The  ICG  measurements  in  equation  (14)  are  collected  at  discrete  time  instances,  t  =  kT,  k  =  0, 1, ..., 
where  T  is  the  sampling  period.  Therefore,  the  continuous  model  described  in  equations  (14)  has  to  be 
discretized.  To  simplify  our  notation,  we  shall  use  C (k)  —  C (kT)  and  m (k)  =  m (kT).  The  discrete  KF 
system  and  the  measurement  models  are  given  as  follows: 


C(A;  +  l)  =  Kd(a„)C(fc)  +  u>(A:) 
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m(k)  =  Vd(oLn)C(k)  +  r)(k ) 


(16) 


where  Kd{ctn)  —  is  the  discrete-time  KF  system  matrix  and  Vd(otn)  —  V(otn)  is  the  discrete-time 

KF  measurement  matrix.  u>(k)  and  rj(k)  are  zero  mean  Gaussian  white  noise  processes  with  covariances 
matrix  Qd  and  variance  respectively.  Discretization  of  state-space  models  can  be  found  in  various 
system  theory  books,  see  for  example  [25]. 

An  explicit  form  of  the  discrete  KF  model  for  the  two-compartment  case  is  given  as  follows: 

Ce(k+  1)  Til  m  Ce(k) 

=  +  io(k)  (17) 

Cp(k  + 1)  r2 1  r22  Cp(k) 


where  is  the  ith  row  and  jth  column  entry  of  the  system  matrix  Kd{a.f).  Note  that  the  matrix  entry 
Tij  is  an  exponential  function  of  the  parameters  k^\  k^  and  k^t. 

To  simplify  the  estimation  process,  we  shall  first  estimate  the  matrix  entries,  r^,  of  the  discrete-time 
system  matrix  Kd(atn )  and  then  compute  the  pharmacokinetic  parameters  for  each  compartmental  model. 

B.  Modeling  of  ICG  Pharmacokinetic  Parameters  and  Concentrations  in  Extended  Kalman  Filter  Frame¬ 
work 

Kalman  filter  provides  a  recursive  method  to  estimate  the  states  in  state-space  models,  in  which  the 
states  are  driven  by  noise,  and  the  measurements  are  collected  in  the  presence  of  measurement  noise 
[26],  [27],  [28].  In  the  case  of  non-linear  state-space  models,  the  extended  Kalman  filter  linearizes  the 
model  around  the  current  state  estimate,  and  then  applies  the  KF  to  the  resulting  linear  model.  The  EKF 
framework  is  also  utilized  for  the  joint  estimation  of  the  unknown  system  and/or  measurement  parameters 
and  states.  In  a  linear  state-space  model  when  both  states  and  system  parameters  are  unknown,  the  linear 
state-space  model  can  be  regarded  as  a  non-linear  model  in  which  the  linear  system  parameters  and 
states  are  combined  to  form  the  new  states  of  the  non-linear  model.  This  system  is  then  linearized  and 
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solved  for  the  unknown  states  using  the  KF  estimator  [29],  [30],  [31].  In  our  problem,  our  objective  is 
to  simultaneously  estimate  the  states,  i.e.,  the  ICG  concentrations  in  each  compartment,  and  the  system 
and  measurement  parameters,  i.e.,  the  pharmacokinetic  parameters  and  the  volume  fractions. 

We  consider  a  Taylor  series  approximation  to  the  non-linear  system  function  at  the  previous  state 
estimates,  and  that  of  the  measurement  function  at  the  corresponding  predicted  position.  This  approach 
provides  a  simple  and  efficient  method  to  handle  the  non-linearity  in  the  new  system  and  measurement 
models. 

Let  6n(k)  denote  the  discrete-time  parameter  vector  of  the  pharmacokinetic  rates  and  volume  fractions 

at  time  k.  For  example,  in  the  two-compartment  model,  02  (k)  is  given  as 

r  i  T 


e2(k) 


(2)  (2) 

Til  T12  T2i  r22  Vl  VKp 


(18) 


In  the  EKF  framework,  Qn{k)  is  treated  as  a  random  process  with  the  following  model: 


en(k  +  1)  =  dn(k)  +  q(k), 


(19) 


where  <;(k)  is  a  zero  mean  white  noise  process  with  covariance  matrix  Sd. 

We  append  the  parameter  vector  On(k  +  1)  to  the  ICG  concentration  vector  C (k  +  1)  to  form  the  new 
non-linear  state-space  model  given  as 


C(fc  +  1) 

K(On)C(k) 

u>(ft) 

(20) 

= 

+ 

6n(k  + 1) 

0n(k) 

c(k) 

m(fc)  = 


Vd(0n)  0 


C(k) 

0n(k) 


+  V(k), 


where  K(0n)  =  Kd(an). 

The  choice  of  Qd,  Sd  and  is  crucial  to  the  performance  EKF  estimator.  It  was  shown  that  if  these 
values  are  selected  less  than  the  actual  values,  it  leads  to  overconfidence  in  the  accuracy  of  the  estimates 
of  the  error  covariance  matrix.  Therefore,  these  matrices  are  regarded  as  tuning  parameters  and  not  as 
the  estimates  of  the  true  covariance  matrices  [32]. 
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C.  EKF  Joint  Estimation  for  ICG  Concentrations,  Pharmacokinetic  Parameters,  and  Volume  Fractions 


In  this  section  we  will  summarize  the  major  steps  of  the  EKF  estimator  for  the  joint  estimation  of 
ICG  concentrations  and  compartmental  model  parameters. 

Let  the  subscript  k\t  denote  the  estimate  at  time  k  given  all  the  measurements  up  to  time  t.  Then  the 
1-step  ahead  prediction  of  the  ICG  concentrations  and  the  compartmental  model  parameters  are  given  as 
follows:  _ 


c 

K(0n)C 

A 

0n 

fc|fc-l 

On 

(21) 


For  the  two-compartment  model,  Equation  (21)  becomes 


Ce 

rjlCe  +  Tl2Cp 

Cp 

— 

T2iCe  +  rhCp 

@2 

fc|fe— 1 

02 

(22) 


fc— i|fc— l 

The  error  covariance  matrix,  Pk |fc_x,  of  the  1-step  ahead  predictions  is  given  as  follows: 


Pk\k-1  ~  Jk-\Pk-l\k-\Jk-\  + 


Qd  0 

o  Sd 

where  Jk  is  the  Jacobian  of  the  non-linear  EKF  system  function  at  time  k.  Explicitly,  it  is  given  as: 


(23) 


Jk  = 


K(0n)  J;[K(0n)C] 


(24) 


k\k 


0  I 

where  0  and  I  denote  zero  and  identity  matrices,  respectively.  The  Jacobian  matrix  for  the  two-compartment 
model  becomes 


Jk  — 


/  „  „  \  / 

Til  T12 

T21  T22 

0(6x2) 


Ce  Cp  0  0  0  0 


0  0  Ce  Cp  0  0 


■■(6x6) 


(25) 


J  k\k 
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where  0(6x2)  is  a  6  x  2  zero  matrix,  and  I(6  x  6)  is  a  6  x  6  identity  matrix. 

The  1-step  ahead  predictions  are  updated  to  the  kth- step  estimates  by  means  of  the  Kalman  gain  matrix 
which  is  given  by 

Gk  =  Pfe|fe-iAT[APfelfc_iAT  +  al}~\  (26) 


where  A  is  the  following  vector 


Vd{0)  ^[Vd(0)C] 
For  the  two-compartment  model  A  vector  becomes 


Jfc|fc-i 


ve(2)  Up(2)  0  0  0  0  Ce  Cp 


k\k-l 


(27) 


(28) 


The  kth-step  estimate  of  the  concentrations  and  the  parameters  are  obtained  recursively  using 


c 

c 

0 

k\k 

0 

+  Gk(m{k)-[Vd{d)C]k\k-1). 


(29) 


For  the  two-compartment  case,  the  kth- step  estimate  of  the  concentrations  and  the  parameters  is 


Ce 

Ce 

Cp 

= 

Cp 

+  Gk(m(k)  -  (ve^Ce  -  vp^Cp)k\k-x)). 

(30) 

0(2)  _ 

k\k 

02 

k\k-l 

The  error  covariance  matrix,  Pk\k>  of  the  kth- step  estimates  is  updated  as 

Pk]k  =  [I-GkA\Pklk_1,  (31) 

where  I  is  the  identity  matrix. 

The  initialization  of  the  ICG  concentrations,  pharmacokinetic  parameters,  and  the  volume  fractions 
plays  an  important  role  in  the  performance  of  the  EKF  algorithm.  Theoretically,  the  state  estimates  can 
be  initialized  at  the  expected  value  of  the  ICG  concentrations,  i.e.  J5[C(0)]. 

One  approach  to  the  initialization  of  the  parameters  is  to  utilize  the  state-space  presentation  given  in 
equation  (16).  Since  E( m(0))  =  Vd(0n(O))-E[C(O)],  m(0)  -  Fd(0n(O))P[C(O)]  is  a  zero  mean  random 
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i 
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variable.  If  we  express  the  variance  of  the  measurement  m(0)  in  terms  of  the  variance  of  C(0)  using  the 
measurement  model  in  equation  (16),  and  solve  for  6n,  we  get  the  estimate  0n(0)  as  the  most  appropriate 
value  for  initialization.  The  details  of  the  selection  of  the  initial  values  for  the  parameters  can  be  found 
in  [27], 

The  initialization  of  the  error  covariance  matrix  is  also  important  for  the  performance  of  the  EKF.  The 
error  covariance  matrix  is  the  matrix  which  provides  information  about  the  error  bounds  for  the  estimates. 
Theoretically,  the  initial  error  covariance  matrix  is  a  diagonal  matrix  where  the  diagonal  entries  are  the 
initial  estimates  of  the  variance  of  concentrations  and  pharmacokinetic  parameters,  i.e. 

Po\o  — 

D.  Compartmental  Model  Order  Selection 

We  adopted  Bayesian  information  criteria  (BIC)  for  the  optimal  model  order  selection.  BIC  is  a  well 
known  information  theoretic  criteria,  in  which  the  optimal  model  order  is  selected  by  minimizing  a 
cost  function  to  avoid  overfitting.  The  cost  function  depends  on  the  number  of  observations,  number  of 
unknown  parameters  to  be  estimated  and  the  maximum  likelihood  function.  A  detailed  discussion  on  BIC 
can  be  found  in  [33],  [34],  [35], 

In  order  to  calculate  the  BIC  for  different  compartmental  models,  we  first  derived  a  likelihood  function 
for  the  extended  Kalman  filter.  The  derivation  is  based  on  the  maximum  likelihood  estimation  of  the 
parameters  in  the  Kalman  filtering  framework  given  as  in  [36],  [37].  We  then  modified  this  likelihood 
function  for  the  extended  Kalman  filter  estimator  for  the  joint  estimation  of  compartmental  model 
parameters  and  concentrations.  The  details  of  the  derivation  is  provided  in  the  appendix. 

IV.  Experimental  Results  -  ICG  Pharmacokinetics  in  Fischer  Rat  Data 

We  applied  the  proposed  EKF  framework  to  the  pharmacokinetic  analysis  of  ICG  data  obtained  from 
four  Fischer  rats  with  adenocarcinoma.  R3230ac  adenocarcinoma  cells  were  injected  below  the  skin  into 
four  Fischer  rats  3  weeks  prior  to  measurements.  The  tumor  size  for  the  rats  ranges  in  diameter  from  5 


Cov{C{  0))  0 

0 


(32) 
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to  30  mm.  The  ICG  concentration  data  was  collected  with  an  MRI-NIR  imager.  The  configuration  of  the 
apparatus,  the  data  collection  procedure,  and  the  details  of  the  experimental  approach  have  been  reported 
in  [9],  [38],  [39], 

Figure  4  presents  the  ICG  concentrations  (/ iM )  from  four  different  rats.  Tumors  in  Rat  1  and  2 
are  classified  as  necrotic  because  of  their  low  tissue  oxy-hemoglobin,  low  total  hemoglobin,  and  low 
gadolinium-diethylene-triamine  penta-acetic  acid  (Gd-DTPA)  enhancement  levels.  Tumors  in  Rat  3  and 
4  are  classified  as  edematous  due  to  their  high  water  content  [40],  It  can  be  observed  from  Figure  4 
that  the  necrotic  cases  display  low  peak  ICG  concentration  values  and  slowly  rising  slopes  unlike  the 
edematous  cases  with  high  peak  values  and  sharp  rising  slopes. 

We  estimated  the  pharmacokinetic  rates  for  the  four-,  three-  and  two-compartment  models.  The  results 
are  given  in  Tables  I,  II,  and  III,  respectively.  The  error  bounds  on  the  estimates  are  derived  from 
the  covariance  matrix  of  the  EKF  estimator.  The  estimated  pharmacokinetic  rates  for  all  compartmental 
models  indicate  that  the  exchange  rates  between  the  capillary  and  the  adjacent  compartment  (ISS  or  EES), 
ka,k%,  n  —  2, 3, 4,  are  significantly  different  for  the  necrotic  and  edematous  tissue.  We  observe  that  for 
the  four-  and  three-compartment  models,  the  estimated  exchange  rates  between  the  ISS  and  parenchymal 
cell  compartments,  A;™,  A$,  n  =  3, 4,  are  comparable.  Similarly,  the  estimated  rate  of  drainage  out  of  the 
plasma,  k™ut,  n  =  2,3, 4,  are  consistent  for  all  models. 

Based  on  the  model  parameter  estimates,  we  computed  the  BIC  values  for  each  rat  data  to  reveal 
overfitting.  The  BIC  values  and  the  number  of  unknown  parameters  for  each  rat  data  are  tabulated  in 
Table  IV.  The  BIC  suggests  that  the  two-compartment  model  is  sufficient  for  all  four  measurement  sets. 

We  further  analyze  the  goodness-of-fit  of  the  compartmental  models  by  means  of  the  residual  analysis. 
The  basic  idea  of  the  residual  analysis  is  to  compare  the  actual  measurements  m (k)  with  its  1-step  ahead 
predictions,  m(A:)fe|fc_1,  based  on  the  estimated  parameters.  A  detailed  discussion  on  residual  analysis 
can  be  found  in  [26],  [41],  The  mean  and  variance  of  the  residual  error  for  the  four-,  three-  and  two- 
compartmental  models  are  tabulated  in  Table  V.  To  normalize  the  error  with  respect  to  the  magnitude 
of  the  actual  measurements,  we  calculated  the  signal-to-noise  ratio  (SNR)  using  the  median  value  of  the 
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measurements  and  the  mean  of  the  residual  errors  for  each  rat  data.  As  seen  from  the  results  in  Table  VI, 
the  SNR  values  are  higher  for  the  two-compartment  case  for  all  data  sets.  These  results  show  that  the 
two-compartment  model  provides  the  minimum  bias  and  the  best  statistical  efficiency.  Figure  5  shows  the 
measured  total  concentration  data  and  its  1-step  ahead  prediction  based  on  the  two-compartment  model 
for  each  rat  data.  Clearly,  there  is  a  good  agreement  between  the  actual  and  the  estimated  measurements. 

Based  on  the  BIC  and  residual  analysis,  we  conclude  that  the  two-compartment  model  provides  the 
best  statistical  fit  for  the  rat  data  and  investigate  the  estimated  model  parameters  in  more  detail. 

In  the  two-compartment  model,  the  rate  of  leakage  into  the  EES  from  the  capillary,  ka  ,  range  from 
0.0106  to  0.0777  sec-1  and  the  rate  of  drainage  out  of  the  EES  and  into  the  capillary,  k[2\  range  from 
0.0247  to  0.0840  sec-1.  Note  that  the  permeability  rates  for  the  necrotic  cases  are  lower  than  the  ones 
observed  for  the  edematous  cases.  Additionally,  the  estimated  values  for  the  pharmacokinetic  rates  are 
much  higher  than  the  normal  tissue  values  due  to  the  increased  leakiness  of  the  blood  vessels  around  the 
tumor  region  [9],  [42].  The  estimated  plasma  volume  fractions  agrees  with  the  values  reported  earlier  [9], 
and  the  values  presented  in  the  literature  [43],  [44].  These  results  confirm  that  v^P  can  be  significantly 
large  in  tumors  and  that  its  magnitude  varies  with  respect  to  the  stage  of  the  tumor  [24].  The  estimated 
values  of  EES  volume  fraction,  v[P  ,  range  from  0.218  to  0.53,  in  agreement  with  0.2  to  0.5  range  reported 
earlier  [23].  Note  that  these  results  are  valid  only  for  the  ICG  pharmacokinetics  in  tumor  cells  R3230ac, 
adenocarcinoma  and  may  not  be  generalized  for  other  types  of  contrast  agents  or  tumor  types. 

Figure  6  shows  the  estimated  ICG  concentrations  in  the  plasma  and  the  EES  compartments  for  the 
two-compartment  model  for  Rats  1  to  4.  Note  that  initial  estimates  of  concentrations  are  noisy  due  to  the 
limited  data  used  in  the  recursive  EKF  estimation.  This  can  be  improved  by  Kalman  backward  smoothing 
[45].  The  peak  values  of  the  plasma  concentration,  Cp,  range  from  2.72  fj,M  to  4.28  //M.  The  absolute 
value  of  the  concentrations  may  not  be  very  useful.  However,  concentration  of  ICG  in  a  compartment 
relative  to  the  one  in  another  compartment  may  provide  useful  information.  We  consider  the  ratio  of  the 
peak  concentrations  in  the  plasma  and  EES  as  a  potential  parameter  to  discriminate  different  tumors. 
The  peak  Cp/Ce  ratio  for  Rats  1  to  4  is  0.551,  0.593,  0.787,  1.151,  respectively.  This  ratio  is  higher  in 
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edematous  cases  consistent  with  the  fact  that  ICG-albumin  leaks  more  into  the  EES  in  edematous  tumors. 


Additionally,  the  ICG  concentration  in  plasma  decays  faster  than  the  ICG  concentration  in  EES  due  to 
its  elimination  through  the  liver  and  kidneys. 

V.  Conclusion 

In  this  paper  we  present  three  different  compartmental  models,  an  extended  Kalman  filtering  framework 
for  the  modeling,  and  estimation  of  ICG  pharmacokinetics  in  cancerous  tumors  based  on  NIR  measure¬ 
ments.  Additionally,  we  introduce  information  theoretic  criteria  and  residual  analysis  for  model  selection 
and  statistical  validation.  Proposed  compartmental  models  are  fit  to  data  obtained  from  Fischer  rats  with 
adenocarcinoma  cells.  The  pharmacokinetic  rates  and  volume  fractions  are  estimated  for  all  models.  The 
estimated  rates  for  all  compartmental  models  indicate  that  the  exchange  rates  between  the  capillary  and 
the  adjacent  compartment  (ISS  or  EES)  are  significantly  larger  for  the  edematous  tissue  as  compared 
to  the  necrotic  cases.  Based  on  the  BIC  and  residual  analysis,  we  concluded  that  the  two-compartment 
model  provides  the  best  statistical  fit  for  the  rat  data  and  ICG  pharmacokinetics.  Parameters  of  this 
model  indicate  that  the  permeability  rates  are  higher  for  edematous  cases  as  compared  to  the  necrotic 
tumors.  Additionally,  we  estimated  the  ICG  concentrations  in  different  compartments.  The  concentrations 
in  different  compartments  may  provide  additional  parameters  for  tissue  characterization. 

While  our  study  indicates  that  two-compartment  provides  the  best  fit  for  the  ICG  pharmacokinetics, 
the  four  or  three-compartment  models  may  be  advantageous  for  modeling  the  pharmacokinetics  of 
functionalized  optical  contrast  agents  that  actively  accumulate  or  activate  in  diseased  tissue  [46],  [47], 
[48].  In  the  near  future,  we  plan  to  investigate  three  and  four-compartmental  models  in  the  EKF  framework 
for  these  optical  agents  collected  from  animals,  and  the  ICG  data  collected  from  human  subjects. 
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Appendix 


The  cost  function  for  BIC  is  given  as 


< Pbic(p )  =plniV  -  21nL(0p,m(l))m(2)) . ,  m(AT)),  (33) 


where  p  is  the  dimension  of  0P,  which  is  related  to  the  number  of  compartments  in  the  model,  N  is  the 

data  length,  and  L(0,  m(l),  m(2), . ,  m (N))  is  the  likelihood  function. 

The  likelihood  function  for  the  EKF  is  given  as 


1  N  i  N 

L(d,  m(l),  m(2), . ,  m  (N))  =  ln[det(fffc)]  - ATkH^Ak, 


(34) 


fe= l 


fc=i 


where  the  matrix  H  is  defined  as: 


Hk  =  APk\k_1AT  +  al, 


(35) 


and  A,  and  Pk\k-\  are  as  defined  in  Section  III.C. 

The  matrix  A  is  defined  as: 

Ak  =  m(k)-[Vd(0)C}k\k_1,  (36) 

where  m(k)  is  the  ICG  concentration  data  collected  from  Fisher  rats  at  time  k,  and  [Vd(9)C]k\k_i  is  the 
1-step  ahead  estimate  of  the  volume  fractions  and  concentrations. 

The  explicit  form  of  the  likelihood  function  for  BIC  calculation  is  given  as 

1  N 

L(0,  m(l),m(2), . ,  m(JV))  =  ^  +  a|)] 

1  k=  1 

H[m(fc)  -  [Vd(9)C]k\k-i]T [APk\k-iAT  +  alr^mik)  -  [V'd(0)C]fc|A;_1] . 

1  fe=i 

where  all  the  parameters  and  matrices  are  as  defined  in  Section  III.C. 
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TABLE  I 


Four-Compartment  Model:Estimated  pharmacokinetic  parameters  using  EKF  algorithm 


*.(-») 

a 

^4) 

jl(4) 

n>c 

*2° 

rve 

L  (4) 

Kj 

jl(4) 

^ out 

(sec-110-2) 

(sec-1 10-2) 

(sec-1 10-2) 

(sec-1 10-2) 

(sec-1 10-2) 

(sec-110-2) 

(sec-1 10-3) 

Rat  1  (Necrotic) 

1.45±0.013 

1.22±0.019 

1.86±0.017 

2.02±0.026 

2.7430.041 

2.41±0.051 

4.05±0.059 

Rat  2  (Necrotic) 

3.4830.048 

2.77±0.034 

4.2830.048 

4.3330.040 

2.98±0.048 

3.03±0.061 

4.76±0.062 

Rat  3  (Edematous) 

4.94±0.052 

5.1630.067 

4.2230.052 

4.13±0.067 

4.14±0.070 

4.27±0.078 

5.3930.085 

Rat  4  (Edematous) 

5.25±0.053 

5.3130.063 

5.0730.068 

5.2230.063 

4.4330.075 

4.03±0.072 

3.853=0.056 

TABLE  n 

Three-compartment  Model:Estimated  pharmacokinetic  parameters  using  EKF  algorithm 


uO) 

/XlQ, 

ju(3) 

Ac 

k? 

1.(3) 

"'out 

(sec-110-2) 

(sec-110-2) 

(sec-1 10-2) 

(sec-110-2) 

(sec-1 10-3) 

Rat  1  (Necrotic) 

1.9330.061 

1.283:0.049 

1.82±0.032 

2.02±0.041 

3.893:0.052 

Rat  2  (Necrotic) 

4.41±0.074 

2.483=0.067 

4.87±0.066 

5.033=0.057 

5.453=0.071 

Rat  3  (Edematous) 

4.71±0.085 

3.883=0.077 

4.953=0.059 

4.68±0.050 

4.423=0.040 

Rat  4  (Edematous) 

5.29±0.091 

6.483=0.096 

4.483=0.062 

4.203=0.048 

5.013=0.055 

TABLE  III 

compartment  Model:Estimated  pharmacokinetic  parameters  and  volume  fractions  usin 

ALGORITHM 

1.(2) 
ib  a 

JU(2) 

Kout 

(sec-1 10-2) 

(sec-1 10-2) 

(sec-1 10-3) 

(10-2) 

(10-2) 

Rat  1  (Necrotic) 

2.47±0.043 

1.063=0.052 

4.613=0.073 

21.83=1.92 

1.413=0.053 

Rat  2  (Necrotic) 

3.543=0.082 

2.983=0.086 

4.833=0.092 

25.43:3.49 

2.423=0.088 

Rat  3  (Edematous) 

6.903=0.101 

4.933=0.072 

3.953=0.048 

30.43=2.81 

4.8430.120 

Rat  4  (Edematous) 

8.40±0.114 

7.773=0.091 

4.023=0.068 

53.03=4.73 

7.0330.321 
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TABLE  IV 


Test  for  the  model  order  selection  for  three  different  compartmental  models  for  four  different  data 

SETS 


Rati 

Rat2 

Rat3 

Rat4 

Model 

P 

4>BIC(P ) 

<j>BIc(p) 

4>BIC(P ) 

<I>BIC(P) 

Two-compartment  Model 

7 

-178.242 

-198.367 

-202.81 

-172.098 

Three-compartment  Model 

11 

-71.615 

-83.849 

-92.182 

-63.912 

Four-compartment  Model 

15 

-39.719 

-45.121 

-56.340 

-30.023 

TABLE  V 

The  mean  and  variance  of  the  error  between  the  estimates  and  measurements 


Four-compartment  Three-compartment  Two-compartment 


Mean 

Variance 

Mean 

Variance 

Mean 

Variance 

Rati 

0.0987 

7.641e-004 

0.0605 

4.732e-004 

0.0072 

2.567e-005 

Rat2 

0.1043 

9.152e-004 

0.0767 

3.017e-004 

0.0057 

4.829e-005 

Rat3 

0.1204 

8.905e-004 

0.0883 

4.921  e-004 

0.0041 

3.021e-005 

Rat4 

0.0904 

5.977e-004 

0.0589 

6.839e-004 

0.0076 

8.618e-005 

TABLE  VI 

SNR  VALUES  FOR  THREE  DIFFERENT  COMPARTMENTAL  MODELS  FOR  FOUR  DIFFERENT  DATA  SETS 


Rati 

Rat2 

Rat3 

Rat4 

Model 

SNR  (dB) 

SNR  (dB) 

SNR  (dB) 

SNR  (dB) 

Two-compartment  Model 

73.2 

68.1 

108.3 

107.9 

Three-compartment  Model 

30.7 

36.1 

23.9 

47.0 

Four-compartment  Model 

20.8 

29.9 

27.7 

18.4 
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Fig.  1.  A  simple  illustration  of  the  capillary  extracapillary  structure. 

Fig.  2.  An  illustration  of  the  ICG  flow  (a)  in  tight  capillary  of  normal  vessel,  (b)  in  permeable  capillary 
of  tumor  tissue. 

Fig.  3.  Block  diagram  of  (a)  the  four-compartment,  (b)  the  three-compartment,  and  (c)  the  two-compartment 
models  for  ICG  pharmacokinetics. 

Fig.  4.  ICG  concentrations  measured  in  tissue  for  four  different  rats. 

Fig.  5.  ICG  concentration  measurement  data  and  1-step  prediction  of  the  measurements  for  four  different 
rats. 

Fig.  6.  ICG  concentrations  in  plasma,  Cp(t)  and  EES,  Ce(t),  for  four  different  rats,  (a)  Rati,  (b)  Rat2, 
(c)  Rat3,  and  (d)  Rat4. 


